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SUMMARY 


The analysis of asymmetric, curved (Coanda) ejector flow has been 
completed using a finite difference technique and a quasi-orthogonal streamline 
coordinate system. The boundary= layer- type jet mixing analysis accounts for 
the effect of streamline curvature in pressure gradients normal to the stream- 
lines and on eddy viscosities. The analysis assured perfect gases, free of 
pressure discontinuities and flow separation. The analysis treats the three 
compound flows of supersonic and subsonic streams, those are: (1) primary 

flow of the driving nozzle, (2) secondary flow between the primary nozzle and 
the Coanda surface, (3) tertiary flow between the primary nozzle and the other 
surface of the mixing section. 

A test program was completed to measure flow parameters and ejector 
performance in a vented Coanda flow geometry for the verification of the com- 
puter analysis. A primary converging nozzle with a discharge geometry of 

3 

0.003175 m x 0.2032 m was supplied with 0.283 m /sec of air at about 241.3 kPa 
absolute stagnation pressure and 82®C stagnation temperature. 

One mixing section geometry was used with a 0.127 m constant radius 
Coanda surface. Eight tests were run at spacings between the Coanda surface 
and primary nozzle 0.01915 m and 0.318 m and at three angles of Coanda turning: 
22.5% 45.0°, and 75.0°. 

The wall static pressures, the locii of maximum stagnation pressures, 
and the stagnation pressure profiles agree well between analytical and experi- 
mental results. 

Key Words : 


Ejector Computer Program 

Coanda Experimental 

Compressible Flow 
Analysis 

Finite Difference 
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Section 1 


INTRODUCTION 

1.1 Background 

The augmentor wing concept under investigation by NASA for STOL air- 
craft lift augmentation Is powered by an air to air ejector. The wing boundary 
layer is drawn into the deflected double flap augmentor channel at the trailing 
edge of the wing and is pressurized by a high velocity slot jet which is ori- 
ented at an angle to the augmentor channel. To predict the performance and to 
optimize the design of the complete augmentor wing, an analytical method is 
needed to predict the performance of the air ejector which powers the augmentor 
flap section. 

Under contract NAS2-5845 a computer analysis was developed for 
single nozzle axisymmetric ejectors with variable area mixing sections using 
integral techniques. Reference 1. The ejectors of primary interest in that 
program and earlier programs were high-entrainment devices using small amounts 
of supersonic primary flow to pump large amounts of low-pressure secondary flow. 
Good agreement was achieved between analytical and experimental results. 

The integral analytical techniques used to analyze the axisymmetric 
ejector configurations are also valid for the analysis of two-dimensional ejec- 
tors. However, the augmentor wing configuration may include asymmetric geo- 
metries, inlet flow distortions, wall slots, and primary nozzles that are at 
large angles to the axis of the augmentor mixing section. The integral tech- 
niques are not easily adaptable to these more complex flows. Finite difference 
techniques can be used to analyze these more complex flow geometries at the 
expense of increased computer time. 

Under contract NAS2-6660 a computer program (Reference 2) was 
developed for two-dimensional, symmetrical mixing sections using finite dif- 
ference technique and rectangular coordinate system. VJhen Coanda effect is 
used in two-dimensional ejectors the geometry is not symmetrical and the use 
of rectilinear coordinates becomes difficult. Flow computations have to account 
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for the pressure gradient In the direction nonnal to the streamlines. 

1.2 Ob.jectlves of the Program 

The following objectives were defined for this investigation: 

1. Develop a computer program for two-dimensional vented Coanda 
ejectors with non-symmetrlc variable area mixing sections, with 
variable Coanda turning, and with variable primary nozzle 
spacing. 

2. Obtain test results with a variable nozzle position vented 
Coanda ejector configuration for the development and checking 
of the computer program. 
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Section 2 


NOMENCLATIURE 
2 2 

Nozzle discharge area (m ; in ) 

Coefficient' appearing in the finite difference equations (-) 

Coefficient: appearing in the finite difference equations (-) 

Specific heat at constant pressure (J/kg*K; Btu/lbm®R) 
Nozzle discharge coefficient (-) 

Coefficient appearing in the finite difference equations (-) 

Coefficient appearing in the finite difference equations (-) 

Nozzle exit height (m; ft) 

Dimensionless eddy viscosity, v^/u^d (-) 

Thermal conductivity of fluid (W/mK; Btu/hrft^F) 

Curvature of streamline, 1/R (1/m; 1/ft) 

Dimensional constant, (32.2 Ibm-f t/lbf-sec ) 

Prandtl mixing length (m; ft) 

Dimensionless mixing length, il/ d (~) 

Node points along a streamline (“) 

Streamline coordinate (normal to streamlines) (m; ft) 
Streamline designation (-) 

Barometric pressure (kPa; psia; inch H 2 O) 



Static pressure (kPa; psig; inch H^O) 

Reference pressure, primary stagnation pressure (kPa; psia) 

Turbulent Prandtl number, v /e (-) 

T n 

Prandtl number, uC /k (-) 

P 

2 2 

Effective heat transfer (between streamlines) (W/m °C; Btu/ft sec®F) 
Radius of curvature (streamline; wall) (1/m; 1/ft; 1/in) 

Gas constant (Nm/kg®K; lbf-ft/lbm®R) 

Richardson number (2u/R) /(3u/3n) (~) 

Streamwise coordinate (along streamlines) (m; in) 

Nozzle coanda wall spacing (m; inch) 

Atmospheric temperature (®C; “F) 

Fluid temperature (°K; °R) 

Maximum fluid temperature at an x = constant cross section (°K; °R) 

Reference temperature, primary stagnation temperature (°K; ®R) 

Velocity in s direction (m/s; ft/sec) 

Reference velocity, u = /R T__ (m/s; ft/sec) 

o g 01 

Unknown velocity at the nth grid point (m/s; ft/sec) 

Velocity in x direction (m/s; ft/sec) 

Centerline velocity (m/s; ft/sec) 

Maximum fluid velocity (m/s; ft/sec) 
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u 

sec 

W 


m 

W 


n 

W 


s 


W 


t 


X 


X 


AX; dX 


y 

Y 


Ay 

3 



y 


”^eff 


T 

W 


Secondary velocity (m/sec; ft/sec) 

Mixing section total flow (kg/sec m; Ibm/sec in) 

Nozzle flow rate (kg/sec-m; Ibm/sec in) 

Secondary flow rate (kg/ sec-m; Ibm/sec in) 

Tertiary flow rate (kg/sec-m; Ibm/sec in) 

Space coordinate in the axial direction (m; in) 

Dimensionless space coordinate in the axial direction, x/d (-) 
Step size in x- direction (“) 

Space coordinate perpendicular to axial direction (m; in) 

Dimensionless space coordinate perpendicular to axial direction, 
y/d (-) 

Dimensionless distance from wall (-) 

Nozzle angle (coanda turning) (degrees) 

Ratio of specific heats, C /C (-) 

P V 

2 

Stream function (kg ra /sec; Ibm/ft sec) 

3 

Fluid density (kg/m ; Ibm/ft ) 

Fluid density evaluated at a reference temperature, and pres 

sure, (kg/m^; Ibm/ft^) 

2 

Dynamic viscosity (Ns/m ; Ibm/ft sec) 

Total effective shear stress (Pa; psi) 

Local wall shear stress (Pa; psi) 




V 

6 

$ 


2 2 2 

Turbulent eddy conductivity, I 3u/9n (m /s; ft /sec) 

2 2 

Eddy coefficient of heat transfer (m /s; ft /sec) 

2 2 

Kinematic viscosity at local temperature (m /s; ft /sec) 

Local wall boundary layer thickness or jet half width (m; in) 

2 3 

Dissipation function (N/m s; Ibm/ft sec ) 
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Section 3 


ANALYSIS OF VENTED COANDA FLOWS IN AUGMENTOR DUCTS 


3.1 Introduction 


This section presents the analysis of asymmetric curved augraentor 
flows which are steady two-dimensional and compressible and for which the duct 
geometry is general. The method extends the work presented an reference (2) but 
a new analysis has been performed and a new program has been developed. The 
analysis employs the finite-difference technique for representing the equations 
of motion for compressible flow. It is essentially a boundary-layer-type jet 
mixing analysis, written in streamline coordinates for ease of computation of 
curvature effects. The effects of streamline curvature on pressure gradients 
normal to streamline and on eddy viscosities are computed. Magnitudes of stream- 
line curvature effects are estimated by using a quasi-orthogonal coordinate sys- 
tem and assumed variation of curvature with distance in the normal coordinate 
direction. 


The analysis treats the mixing of three compressible flows of the 
same perfect gas under the assumption that initial conditions are known and that 
pressure discontinuities and flow separation are absent. The nozzle exit flow 
may be supersonic but it is assumed that expansion or recompression outside the 
nozzle, if needed, will bring the nozzle stream to the local ambient pressure so 
that shocks and expansion waves at the nozzle exit plane are avoided. Previous 
work has shown that augmentor performance is little affected by moderate degrees 
of departure from conditions of correct nozzle expansion. The flows considered 
include compound flows of supersonic and subsonic streams; however no provision 
is made for compound choking which may occur with an appropriate transverse dis- 
tribution of Mach number. Such a condition is amenable to analytical treatment 
under simplified circumstances, but has not been encountered in experimental 
tests carried out so far. 

To retain the simplicity and speed of the boundary layer approach to 
augmentor calculation, while incorporating approximate curvature effects, it is 
necessary to assume an approximate starting line. In the present work the 
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starting line is comprised of two circular~arcs (See Fig. 1 ) which are tangent 
to each other, and perpendicular to the nozzle axis at the nozzle exit plane; 
one arc is normal to the upper wall and one to the lower wall of the duct. In 
the absence of detailed experimental information on velocity profiles in the 
wall boundary layers and in the jet shearing layers at the initial plane, the 
initialization condition has assumed uniform stagnation pressure in the nozzle 
flow and, separately, for the secondary and tertiary flows. In addition, for 
the examples worked out in this report , the secondary and tertiary flows have 
been assumed to have the same stagnation conditions. In computing initial con- 
ditions around the circular arcs, the effect of curvature on normal pressure 
gradient has been taken into account; the initialization satisfies the conti- 
nuity equation separately for primary, secondary, and tertiary streams under 
the constraints of local duct width (along the assumed circular arc starting 
lines), location of the nozzle centre line, and angle of the jet axis with res- 
pect to the coordinate system of the duct walls. This initialization is of 
course approximate but is reasonable to use in the absence of better informa- 
tion on flow starting conditions. If better information is available the ini- 
tialization procedure adopted in this analysis may readily be replaced to use 
more detailed or exact information. 


3.2 Equations of Motion 


The momentum, normal pressure gradient, 
streamwise coordinates are: 


lii = _ + J_ 

3s 9s 9n 


ef f 


) 


P 


2 

u K = 


9P 

9n 



pu 


9(C T) 

P 

9s 


9P , 9 

9s 9n 


(q 


ef f 


) + ^ 


and energy equations in 


( 1 ) 

( 2 ) 

(3) 


in which 


T 


eff 


(y + pv^) 


9u 

9n 


8 





(k + pc e ) 
P « 


9T 

9n 


$ = (y + pv^) 


W 


2 


In these equations s and n measure distance along and normal to the streamlines, 

respectively, and u is the velocity component in the stream direction; p is the 

static pressure, p the density and T the temperature of the fluid, t __ is the 

erf 

total effective shear stress and the eddy viscosity of the fluid with y being 
the dynamic viscosity. Correspondingly, is the effective heat transfer be- 
tween streamlines with e be.. eddy coefficient of heat transfer and k the 

H. 

fluid conductivity. Constant values of laminar and turbulent Prandtl numbers 
have been assumed in the analysis. The term $ is the dissipation function, in- 
cluded in the energy equation. The first order effects of curvature on static 
pressure are included through the normal pressure gradient equation (2) . 


Stream Function 


The stream function 'F is defined by 


9n 


= pu 


(4) 


With (4) equations (1), (2), and (3) becoire: 


in = ^11^ + 

^ 3s p 3s 


U pu(u + pv.j) If] 


uK = 


3P 


u 


3’F 

3(C T) 

2 

3s 


u _3P _3_ 

p 3s ^ 3'F 


3T 

pu(k + CpPPj,) 


r -t 2 


r 9^1 


3.3 Dimensionless Parameters 


Each variable in the equations of motion is normalized by use of 
the following reference variables: 


* 


u 
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* s 


s = -- 


* 


n 


n 


d 


* 

R 


R 

d 


■k 

K 


d 

R 


* 


y 


y 


P0l“o‘* 




P 

r 



* 

'P 





in which 


u = A 
o g 01 


PQI = 


d = nozzle exit height (the small dimension) 
= primary stagnation pressure 
= primary stagnation temperature 


Introducing these dimensionless groups into the equation of motion 
yields the following results: 

k 


* du* 1 9P* . * 9 

9s p 9s 9¥ 


k k k * - 9u 
P u (y + P E) ^ 


* * 3P 
u K = — r 


94' 


u 


* 9T 


9s 



* k 


f k \ k 

Y - 1 

u 9P , * 9 

* * 

y . E 1 9T 

. y J 

* * ^ * 

P u 

P P J . * 

p 9s 94' 


^ r rt^ 94' _ 


/* - \ 

f k * s 


* 

y - 1 

y + p E 


* * 8u 

y J 

* 

^ P ^ 


^ ^ * 
34' _ 


From henceforth we omit the superscript * for convenience so that 
the following are the equations of motion in dimensionless form: 


u 


9s 


LIZ + _L 

p 9s 94' 


P 


pu(y + pE) 


9u 

94' 


(5) 
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( 6 ) 


uK = 


8P 

av 


9T _ 

fr - 1] 

u ap a 

pu 

f y ^ E ] 

3T 

ds 

L T J 

p as “ a'F 


a'i' 


fy ” i1 

V + pE] r 9u 

1 Y J 

P J L 91 


(7) 


3.4 Evaluation of the Eddy Viscosity 


' In general the eddy viscosity is defined by: 
.2 9u 

'■t = ^ 


( 8 ) 


Writing this in dimensionless form, using the dimensionless parameters specified 
previously and the stream function, the eddy viscosity expression becomes: 


E = L pu — 


(9) 


in which 


L 


d 


and E = 



The effect of streamline curvature is taken into account by use of 
the Richardson Number correction in the following approximate form: 


L = L exp(-3R.) 

O 1 

R. > 0 
1 

(10) 

L = - expORpl 

R. < 0 
1 

(11) 


in which is the dimensionless mixing length in the absence of stream curva- 
ture and is the Richardson Number defined by: 


or 


R, 


2u , 3u 
R ' an 

2K 



( 12 ) 
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For small values of the above dependence of L on is approx- 

imately in accord with the linear relationship derived by Bradshaw (3) . For 
large values of |R^| an empirical correlation is not available, and the expo- 
nential relationship has been assumed. 

In the absence of curvature the mixing lengths L are defined as 

o 

follows . 


Boundary Layer 


In the inner part of the layer the Van Driest approximation is used. 


= 0.41 Ay[l - exp (-y"^/26) ] 


(13) 


in which Ay is the dimensionless distance from the wall. The variable 


+ Ay^l^w 


(14) 


is evaluated using 

yu. 


Ay 


w Ay, 


_2 ^ 
2 9s 


(15) 


in which the subscript 2 denotes the streamline coordinate point closest to the 
wall . 


In the outer part of the layer the mixing length is evaluated by 


L = 0.09 
o 



(16) 


in which 6 is the boundary layer 99% thickness and d is the nozzle exit width. 


In the middle part of the layer the smaller of the values of 
provided by Equations (13) and (14) is used. 


Jet Shear Layers 

In the shear layer adjacent to the potential-core zone of the pri- 
mary jet the mixing length is evaluated from 
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(17) 


L 

o 


0.07 




u 


1 + 0.6 


sec 


U 


CL 


in which A is the shear layer width (including the zone between 1% and 99% of 
the total velocity difference between primary and secondary streams) divided by 
the nozzle width at its exit plane. The term in parentheses takes approximate 
account of the effects of the secondary velocity, ^ co~flowing outer 

stream, on the mixing of a jet whose centerline velocity is U . 

dj 


For a "fully-rounded" portion of the jet flowing co-axially with a 
secondary potential stream, the mixing strength has been calculated from 


L 

o 


0.09 


U 


1 + 0.6 


sec 


CL 


(18) 


\ 

in which A is the half-width of the jet (evaluated from centerline to the point 
at which difference between local and secondary velocity is only 1% of the dif- 
ference between centerline and secondary velocity) , divided by the nozzle width 
at its exit plane. j 

Developing Pipe Flow Region , 


For the region downstream at the point where the jet spreads to 
intersect the edge of the boundary layer the mixing is evaluated, as a first 
approximation only, from 



f 

1 — 

. -t2 I- 

1 ^ 

4' 

L =T 

0,14 - 0.08 h 

. - - O.Obfj 

o d 

L 

wj L 




(19) 


in which w is half the total width, and Ay the distance from the wall. This 
formula is due to Nikuradse and is cited by Schlichting (3) for fully developed 
flow in round tubes. Near the wall the mixing length is evaluated by the Van 
Driest approximation cited earlier, provided the local mixing length so calcu- 
lated is less than that given by the Nikuradse formula. 
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The laminar dynamic viscosity is evaluated from (1): 


y 


^ref 

T T 
01 

1/2 

pref + 

198. 7°R 


T 

_ref 



198. 7°R 


( 20 ) 


in which is the laminar dynamic viscosity at a temperature and a pres- 
sure Slid T is the dimensionless temperature (normalized by the inlet sta- 

gnation temperature of the primary stream 


The wall shear stress (normalized by the inlet stagnation pressure 
is given by Equation (15). For this equation to be realistic, the grid 
spacing must be chosen so that y at node point 2 is not greater than 3-4. 
The associated wall friction velocity (normalized by the reference velocity u^) 
is given by: 


^ ii 

Ay 2y dxj P 


( 21 ) 


3.5 Finite Difference Procedure 


By the finite difference technique the derivatives in the differ- 
ential equations of motion are replaced by differences either along a stream- 
line between two neighboring points X and X + dX or normal to it between two 
neighboring points 'I' and '1^ + d'F. The finite difference equivalence of the 
equations of motion are obtained as follows with reference to the following grid 
lines : 
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Writing the velocities in terms of a taylor expansion: 
n n 

— OU A ^1/ -1_ O U 2 

U _ = U ~ At + 

n - 1 n aw 2 2! 

n n 


Eliminating 


2 2 2 9 

^2 AW^^ AW2 AW^^ 

~ % + 1 ~ % - 1 ^ '"n r~ % 


2 2 
AW- AW, 

+ — AW 1- AW — — 

aw 1 2 ^2 2 


Dividing by AW^^AW 2 


AW^ n + 1 AW 


AW AW 

2 % - 1 = Af^ “n - “n + l?|^ f \ + *'' 2 ] 


This equation leads to: 


aw AW, (AW, + AW-) ^ n + 1 
n 11 2 


A 1 ' 2 (A' 1 'j^ + Af^) ^“n “n - 


+ 1 ” + S (u - U ) 

OT Dn + 1 n 4n n-1 


3? = + 1 - - 1> 

n 


in which 


and Sj Af (Af + A¥ ) 
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Using another Taylor Series expansion with a general coefficient 
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Solving for — 


s — 


we obtain: 
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5u 
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(AH'^ + AH'^) 
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(u - u ) 
n n - 1 
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The terms in Equation (5) become: 
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1 
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in which S = pu(y + pE) 
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With these finite-difference equivalents for the derivative terms 
Equation (5) may be written in the form: 


A -u„ +B u u* 

n-12, n n-12, n + 1 n-12, n-1 n-1 


( 22 ) 


in which 


u 

A 1 = .A ^ 
n-1 AS 8 9 

n 


(23) 


®n - 1 = -^8 


^n - 1 ^9 


D 


n-1 


1, n 


3x 


u 


1» n 
m + 1/2 ^ ^n 


(24) 

(25) 

(26) 


u_ (S . ^ + S ) 

, _ 1, n n + 1 n^ 

8 (A’i'^ + A'F^) A'F^ 


(27) 


Ut (S + S J 
y _ 1 , n n n-1 

9 (A*?^ + AH-^) AH '2 


(28) 


S = pu(y + pE) 


(29) 


Similarly, Equation (7) may be written in the form: 


A + B T„ 1 + C -T. 1 = D 

n-12, n n-12, n-1 n-12, n-1 n-1 


(30) 


in which 


u 

A = — =-* — + Y ’ + Y 

n-1 AS 8 ^9 

n 


(31) 


- 1 = -’^8' 


- 1 = -^ 9 ' 


(32) 

(33) 
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D = "l. " T + ^ ” fill 

” - ^ „ l3xj^ ^ 


+ ^ ~ u, S 
Y 1, n n 


®5^“l, n + 1 - "l. 


+ S4< 


"X, n - "l. n - lU 


S’ . T + s' 

y I _ n + 1 n_ 
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^9* 
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JL + -£1l. 
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(34) 

(35) 


(36) 

(37) 


Coefficient Matrix 
The general equation 


A .X + B ,X . , + C ,X = D , 

n~ln n-ln + l n-ln-l n-1 


with boundary conditions 



U = 0 
n 


T = T 
n n - 1 


can be written for each of the n grid points with the result: 
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1 -6 
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2 2 2 


3 


2 




s 
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where 6 = 0 in momentum equation 
6 =» 1 in energy equation 


in which 


The first and second equations are; 


Xi - 6X2 = 0 


"A ^ = “1 


A- 1X2 + 3^X3 = 


A'l = C^6 + 


The last two equations are: 


C „X _ + A 
n - 2 n - 2 n 


nX , + B 
2 n - 1 n 


2 n n - 1 


-6X + X =0 

n - 1 n 


Combining these 


in which 


C „X o + A* 9^ 1 = D 1 

n-2n-2 n-2n-l n-1 


A' « = A - + B _6 
n-2 n-2 n-2 


With these two results the order of the matrix can be reduced to n - 2 . 
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[k\ B, 0 - - 
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This is solved by the Thomas Algorithm as indicated in References 

(2) and (5). 

3. 6 Boundary Conditions 

The boundary conditions of the walls are: 

y = (upper and lower walls) 

K = K (x) (wall curvature) 
w 

Y = const 
u = 0 
9T 

= 0 (adiabatic wall) 

In the program provision is made for the calculating wall curvatures 
from wall coordinates x, y using a least-squares smoothing procedure. Wall co- 
ordinates must be specified with sufficient precision to estimate realistic 
values of wall curvature. 


3. 7 Quasi-Orthogonal Coordinate System 


To preserve orthogonality of the coordinate system the step sizes 
AS must be adjusted according to: 


or 


9As _ A^ 
9n R 

9As _ As 
9W ” puR 


(38) 
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Starting with an arbitrarily chosen step size at the nozzle centre line the 
step size for adjacent streamlines is calculated with the finite-difference 
equivalent of this formula, using mean values of curvature, density, velocity 
and step size between the neighboring streamlines. 

Streamline curvature (K = 1/R) is assumed to decay exponentially 
with distance from the wall according to: 

r 

K = Kj^ll - -^Jexp(-o|K^|toj^) + -j^xp(-a|K2|in2) (39) 

in which (positive) and (negative) are the curvatures of the lower and up- 
per walls, respectively, An^ and An^ are the distances along the streamline 
orthogonal to the lower and upper vjalls, respectively and w is the sum of the 
two. The decay constant a appears to have a value 2 for potential flox^r around 
a cylinder, but in these boundary layer flows a value of 4 seems more appropri- 
ate, from comparisons of measured and calculated pressure distributions. The 
contribution of the first term to streamline curvature in the vicinity of the 
upper wall is negligible, and w. 

3.8 S olution Procedure 


The first step in the solution is to determine all flox\r properties 
on the initial line. First the radii of curvature of the tx\ro circular arcs are 
determined iteratively to satisfy the conditions of mutual tangency at the noz- 
zle centerline, and orthogonality to the nozzle centerline and to the tx*TO walls. 
The initial stagnation conditions are given and assumed the same for secondary 
and tertiary flows (See Fig. 1 ) . The primary mass flow and the sum of primary 
and secondary mass flow rates are also given. The initialization then solves 
iteratively for the split betx^een secondary and tertiary flow rates, and the 
values of all properties along the initial line under the requirement of mass 
conservation, and assuming isentropic flow up to the initial line. The solu- 
tion procedure accommodates substantial effects of streamline curvature, but 
will not succeed if the total ingested flow rate is so small that the local 
velocity e.g. at the lower side of the nozzle (Fig. 1 ) is negative. The 
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initialization procedure also selects the location of streamline node points, 
T^th veiry close spacing required near the walls, and moderately close spacing 
in the initially- thin shear layers of the jet. 

For a set of n streamlines and known boundary conditions, equations 
(22) and (30) each provide a set of n - 2 conditions to solve for the unknovm 
velocities and temperatures. Each set of equations can be solved simultaneously 
if the pressures at dovmstream node points are kno^m or asstjmed. For calcula- 
tion pf flow betvjeen curved channel walls, the pressure gradient along one 
streamline is assumed and the downstream pressure on that streamline determined 
from this gradient and the arbitrary step size. Pressures at corresponding 
node points on adjacent streamlines are determined by use of the finite-dif- 
ference form of Equation (6). With the doxrastream pressures determined, equa- 
tions (22) and (30) are solved to provide downstream velocities and temperatures, 
and subsequently all other properties at downstream .points. If the calculated 
value of the outer boundary location does not agree satisfactorily with the 
actual wall geometry, a new value of the pressure gradient is chosen. 
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Section 4 


TEST PROGRAM 


A two dimensional experimental rig was designed, fabricated, and 
Installed in our laboratory. The purpose of the experimental work was to obtain 
test data for verification and adjustment of the computer analysis. The experi- 
mental program is described in this section. 

4.1 Experimental Apparatus 

4.1.1 Two Dimensional Ejector 

The two-dimensional ejector consists of a slot t 3 rpe primary nozzle 
and a two dimensional mixing section. The arrangement of the ejector system for 
the four nozzle positions tested is shown in Figures 2 and 3 and the positions 
are listed below. 

Test if Nozzle Angle (g) Spacing (t) (inches) 


1 

& 

2 

22.5° 

0.80 

3 

& 

4 

45.0° 

0.80 

5 

Sc 

6 

45.0° 

1.31 

7 

Sc 

8 

67.5° 

0.70 


The nozzle angle (3) is defined as the angle measured between the 
vertical and the line running from the center of the coanda arc to the center- 
line of the nozzle at the throat while the spacing (t) is defined as the per- 
pendicular distance of closest approach between the nozzle and the coanda sur- 
face. 


A picture of the primary nozzle is shown in Figure 4. The discharge 
slot is 0.1215” + .0005" by 8,00" with rounded corners. The side walls are 
made from one quarter inch carbon steel. Four internal supports prevent sub- 
stantial widening of the discharge slot when the nozzle is pressurized. Dial 
indicator measurements performed in previous tests revealed that the slot opened 
up by about 0.0008 inches in the center of the nozzle, about 0,0004" at the 
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quarter X'^idth location, and zero near the ends of the slot. This corresponds 
to an increase in nozzle slot area of 0.33% when pressurized. Stagnation pres- 
sure measurements were made with a kiel probe from side to side at two different 
axial locations and x^ere found to be uniform across the 8" width of the slot 
(See Appendix A) . 

Aluminum pieces on which the nozzle pattern had been cut at the 
desired angle were bolted to the inside of the side plates to position the noz- 
zle in the inlet to the mixing section (Figs. 0, 7, 8). The nozzle was held 
firmly in place by threaded rods connecting brackets on each side of the nozzle 
to brackets attached to busses welded to the side plates of the mixing section 
(See Figs. 6, 7, 8), 

The mixing section, as shoxm in Figure 2, consists of a rectangular 
variable area channel formed by tx*70 identically contoured aluminum plates, tx-70 
flat side plates and two dissimilar curved inlet pieces. The upper bellmouth, 
the coanda surface, is a constant radius (r = 5.00 inch). The nozzle is posi- 
tioned close to this surface and relies upon the coanda effect to turn the pri- 
mary floxj smoothly into the test section. The loxjer inlet consists of a bell- 
mouth piece and a straight section. The pictures in Figures 6, 7, and 8 show 
three views of the mixing section. The two contoured plates were positioned in 
symmetrical locations about the centerline to form the channel tested (throat 
height of 1.875"). The width of the mixing section is a constant 8,00 inches 
along the entire length. The variation of channel height with distance from 
the nozzle discharge is given in Table 1. Three plexiglass windows were in- 
stalled along each side of the mixing section so that tufts of wool mounted 
inside could be observed for indications of flow separations and unsteadiness. 

The screened mixing section inlet is shoxra on Figure 9. Earlier 
tests (Ref. 2) without the extended inlet shox<red that highly sx7irling comer 
vortices were formed in the four comers of the bellmouth and extended into 
the test section. The extended inlet eliminated the comer vortices and im- 
proved the stability of the ejector flow and static pressures. Four sets of 
screened inlets were used in the testing to accommodate the four nozzle positions. 
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4,1.2 Facilities for Ejector Tests 


Three subsystems are required for the operation, control and measure- 
ment of the air flow through the two dimensional ejector (Figure 9). These 
three subsystems referred to as the primary flow, the mixed flow, and the bound- 
ary layer suction systems are discussed below. 

The primary air flow is supplied by a 900 scfm non lubricated screw 
compressor at 100 psig and an equilibrium operating temperature between 100®F 
and 140 ®F. The primary air flow rate and pressure are controlled by an auto- 
matic pressure regulator capable of maintaining pressure to within + .1 psl of 
a set value. The mass flow is measured by a standard 3 inch Denial orifice 
system. A flexible hose connects the primary orifice system to the nozzle. 

The mixed flow system consists of a plenum chamber and an 8 Inch 
orifice system. Two different operating flow rates were achieved by the fol- 
lowing equipment combinations. 

1, Maximum Flow Rate at Atmospheric Discharge - Mixed flow dis- 
charges directly into the laboratory. 

2 . Reduced Flow Rate at Back Pressure - The plenum and orifice 
are connected to the mixing section discharge. 

Mixed orifice flow rates were obtained only for the reduced flow rate condi- 
tions. The plenum is shown connected to the mixing section by flexible hose 
in Figure 10. 

The suction system removes the boundary layer flow from each of 
the four corners of the mixing section to prevent wall boundary layer separa- 
tion in the ejector. Figure 11 shows six 3/4 inch tubes connected to the top 
corners of the mixing section. A total of 12 tubes (top and bottom) collect 
the boundary layer flow from the four corner suction slots which are 0.060 
inches wide and are machined into the sides of the contoured plates (See Figs. 

12 and 13) . The four tubes at one X location are connected to a single large 
tube under the mounting table (Fig. 14). The three large tubes are each 
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connected to a large tank plenum through a separate throttle valve. A Roots 
blower draws the air through the suction system and through a three inch 
orifice system. The suction system is capable of removing about 1% and 2% of 
the mixing section flow rate. During the operation of the ejector rig, the 
boundary layer suction system was used to prevent flow separation in the mixing 
section diffuser. 

The ejector system was operated by starting the primary air flow 
at low pressure and flow rate. The primary nozzle pressure was increased to 
22 psig and the suction was then turned on. Approximately 1/2 hour of warm up 
time was allowed before testing. 

4.2 Instrumentation and Data Reduction 


4.2.1 Instrumentation 


The following instrumentation was used in the test facility. 
Primary Flow System 

Flow Rate ~ Standard 3" orifice system 

Nozzle Pressure - Bourdon Pressure Gage accurate to + .10 psig 
Nozzle Temperature - Copper Constantan thermocouple with digital 

readout 


Mixed Flow System 

Flow Rate - 8" orifice system for reduced flow rate conditions 
(Tests 2, 3, 6 and 7) 

Static Pressures “77 wall static pressure taps located throughout 

the mixing section on both top and bottom con- 
toured plates and on both bellmouth pieces (See 
Figs. 12 and 13). As shown in Figure 10, the 
static taps were connected to a valving system 
which permitted easy determination of individual 
static pressures without the need for a large 
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manometer bank. Tygon tubing was used to con- 
nect the taps to four pressure sampling valves 
each capable of handling 24 inputs. These pres 
sure valves were in turn connected through a 
switching network to any of 3 well type mano- 
meters which permitted the accurate determina- 
tion of pressures over the range of + 25 inches 
of water gage to - 75 inches of water gage. 

Traverse Data - Stagnation pressure and temperature profiles were 
measured at up to 11 axial locations using a 1/8" 
diameter stem kiel pressure- temperature probe. 

Both a mercury manometer and a pressure transducer 
coupled with a direct digital readout were used 
for pressure measurements. A direct digital read- 
out was used to indicate total temperatures . 

Suction Flow System 

Flow Rate - 3" orifice system 

Suction Pressure - Bourdon type pressure gage 

4.2.2 Data Reduction Procedures 


Three types of data reduction calculations were needed in this pro 

gram. 

1. Standard orifice calculations 

2. Velocity profile calculations 

3. Integration of velocity profiles to calculate flow rate 

All of these calculations were programmed on a time sharing com- 
puter. The orifice calculations were programmed as a subroutine to the main 
data reduction program using standard orifice equations and ASME orifice coef- 
ficients. The equations used in the determinations of velocity were included 
in the main program and are the standard compressible flow relationships which 
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can be found In most fluid mechanics text books. 


The integrated mass flow rate for each traverse was computed by 
integrating the product of the local velocity and local density over a two di- 
mensional section of unit width. The program also calculated the "mass-momen- 
txjm” stagnation pressure at each traverse section using the equations presented 
on pages 52 and 53 of Reference 6. The mass momentum method determines the 
flow conditions for a uniform velocity profile which has the same integrated 
values of mass flow rate, momentum, and energy as the non-uniform velocity pro- 
file actually present. 

4.2.3 Experimental Uncertainty 

Orifice Calculations 


The techniques presented in Reference 7 were applied to the primary 
flow orifice calculations and the mixed flow orifice calculations. The follow- 
ing uncertainty results were obtained: 


Orifice 

Primary Nozzle 
Mixed 


Pressure (psig) 
22 psig 

slightly above 
atmospheric 


Uncertainty 
+ 0 . 8 % 

+ 1.3% 


Static Pressure 


Uncertainty in the wall static pressures occur mainly because of 
fluctuations in the manometer liquid columns caused by unsteadiness in the 
flow. The degree of these fluctuations may therefore be used as an indication 
of the uncertainty of the pressure readings. For the unrestricted maximum 
flow rate condition the wall static pressure fluctuation reached a maximum of 
+ 1.0 inch of water, while for the reduced flow rate condition the maximum 
reached only +0.4 inch of water. 

Integrated Mass Flow Rate 

The mass flow rate calculated by integrating the results of the 
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stagnation pressure and temperature traverses is influenced by many items and 
is therefore very difficult to estimate. The following items all contribute 
to the uncertainty in integrated mass flow rate: 


1. unsteady wall static pressures 

2. unsteady traverse stagnation pressures 

3. instrument accuracy of the pressure transducer and digital 
readout 

4. inaccuracies due to the effect of steep velocity gradients on 
sensed pressure 

5. inaccuracies due to probe effect near the mixing section walls 

6. inaccuracy in probe position 

7. assumptions and inaccuracies associated with the data reduc- 
tion computer program 

8. data recording errors or computer data input errors 

9. errors caused by loose connections in the pneumatic sensing 
tube between the probe and the transducer 

10. non-two-dimensional flow distribution across the width of the 
8 inch mixing section 

All of these effects could combine to give both a + uncertainty band and a 
fixed error shift. 

One measure of the uncertainty due to these effects is obtained 
from the limits of individual integrated mass flows for each test run. These 
values are listed on Table 2 for all of the test runs with traverse data. 

The results presented on Table 2 show an average variation of +4.4 and -3.0 
or a total spread of 7.4%. These values only include the effect of variable 
uncertainty and exclude the uncertainty due to probe errors in steep gra- 
dients and near walls and integration assumptions. Both of the excluded er- 
rors probably cause the integrated mass flows to be too large because the 
probe tends to measure too high near the wall and the integration program 
neglects wall boundary layers. > 
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4.3 


Test Schedule and Results 


A total of eight ejector tests were carried out on one mixing 
section configuration (1.875" height) and at one nozzle pressure, 22 psig. 
Four different nozzle positions were tested each at atmospheric discharge and 
at a back pressure condition. Figure 3 summarizes the conditions for each of 
the tests. In each test, readings were taken from the wall static pressure 
taps, orifice system instruments, and from the total pressure and temperature 
taps on the traversing probe. The number and location of traverses varied 
from test to test. The traverse locations and a summary of the data presented 
for each test is given in Table 3. 

The data presented in this report falls into the following cate- 
gories: 

Test Conditions and Mass Flows 

Static Pressures 

Maximum Local Pressures 

Velocity Profiles 

Richardson Number Coefficient Sensitivity 

Streamline Curvature Decay Sensitivity 

A summary of the figures and tables used to present data from each 
test run is presented in Table 3. Discussion of the data will be taken up in 
the next section. 

A sample of the static pressure data as taken is tabulated in Ap- 
pendix A for Tests 7 and 8. 
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Section 5 


COMPARISON OF ANALYTICAL AND TEST RESULTS 

5.1 Test Conditions and Mass Flows 

Table 5 and Figure 15 show a comparison of 

a) the total flow rate measured by an orifice downstream of the 
diffuser 

b) total flow rates determined by integration of measured velocity 
profiles at various sections in the test section 

c) flow rates inferred by use of the computer program with best 
fit to the static pressure in the throat region 

In general the flow rates determined by methods (a) and (c) agreed 
satisfactorily within approximately 6%. Flow rates estimated by integration of 
experimental velocity profiles at various stations were consistent with one 
another within approximately 8% but disagreed with the results of (a) and (c) 
up to 15%. Previous experience (Ref. 2) also showed that integration of experi- 
mental velocity profiles yielded too high a mass flow. In that case the dis- 
crepancy was of the order of 6%. The velocities were determined experimentally 
by the use of a Kiel probe to determine a local stagnation pressure coupled with 
the assumption that the local static pressure was equal to the wall static pres- 
sure. The velocities were determined in this way only for stations downstream 
of high wall curvature. The disagreement between total flows determined by in- 
tegrating velocity profiles and those obtained from orifice measurements may be 
due to the effect of high shear and turbulence level in these flows upon the 
apparent stagnation pressure reading of the Kiel probe. In view of these dis- 
crepancies reliance was placed in these tests upon the mass flows determined 
by methods (a) and (c) . 

The accuracy of the primary flow measurements determined by orifice 
readings for the primary flow are of the order of 3%. The secondary flow de- 
terminations by orifice have an apparent uncertainty level of + or - 1.5%. 
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5.2 


Wall Static Pressures 


Figure 16 shows a comparison between the static pressures predicted 
using this computer program and those measured with the symmetrical diffuser 
employed in the NAS-50 program (Reference 2) for which wall curvatures were 
very small and had little effect upon the axial static pressure profile. This 
comparison which is included as a check point in this discussion shows that 
the curvature program predicts the wall static pressures reasonably well when 
the nozzle is located at the mid-plane of a symmetrical test section. 

Figures 17 through 20 show experimental values of the wall static 
pressures measured for eight experimental cases with the unsymmetrical test 
section and with the nozzle located at various distances from the curved wall 
and at various angles with respect to the test section axis. In general there 
is a considerable difference in static pressure between top and bottom walls 
with the lowest pressure being near the top wall which had the highest cur- 
vature and near which the nozzle was located. These low pressures between 
the top wall and the nozzle were accompanied by velocities considerably higher 
than those in the region between the nozzle and the lower wall. Downstream of 
the region of considerably high wall curvature the measured wall static pres- 
sures were nearly the same on top and bottom walls. 

Figures 17 through 20 also show the computed static pressure dis- 
tributions on the top and bottom walls. In general the agreement between ana- 
lytical measured results is considerably better with high coanda effect, i.e., 
large turning angle (up to 67.5°) and small spacing between the nozzle and 
the wall. The greatest discrepancies between analytical and experimental re- 
sults are associated with those cases where the calculation method indicates 
that the flow is on the verge of separation, for example, cases 1 and 2. In 
cases 5 and 6 the static pressure distribution on the top wall near the noz- 
zle is quite different from the calculated value. Here the spacing between 
nozzle and wall is large at 1.10 inches. This large venting of the flow be- 
tween the nozzle and the wall appears to have substantially diminished the 
Coanda effect lessening the tendency of the flow to cling to the upper wall 
and increasing the possibility of separation in the region immediately 
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downstream of the nozzle. In the region of separation the flow calculation be- 
comes somewhat uncertain and experimental details were insufficient to ascertain 
whether the flow were actually separated in that region. However the degree of 
agreement between measured and calculated results was substantially poorer for 
cases 5 and 6 with large venting between the nozzle and the wall. 

In general the reasons for differences between the analytical and 
calculated results are as follows; 

1. Flow Separation - The experimental results for cases 1 and 2 
and case 6 appear to be on the verge if not actually past the 
margin of flow separation, as indicated by the calculated 
values of the wall shear stress or the velocities near the wall 
for those cases. 

2. Initialization Approximations - As explained in the earlier 
section of the report, the initialization process assumes a 
circular arc starting line for each of the two regions between 
the nozzle and the upper and lower walls. The initialization 
process is assumed to have isentropic flow up to the starting 
linewherethe nozzle has zero thickness, so an approximation 
is used for estimating the decay of streamline curvature with 
distance away from the wall. In the initialization process it 
was found that the calculated wall static pressure distribution 
immediately downstream of the nozzle was very sensitive to the 
ratio of the flows between nozzle and upper and lower walls 
respectively. This flow ratio was not available experimentally 
and was determined in the initialization process by requiring 
smooth continuity of static pressure along the circular arc 
starting lines. The initialization process was also sensitive- 
ly dependent upon decay of wall curvature away from the upper 
and lower walls, especially in the jet zone. 

3. The Use of a Quasi— Orthogonal Coordinate System - In the cal- 
culation method curvatures were estimated by the use of 
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equation 39. In principle this estimation could have been used 
for a first approximation to determine the entire velocity field 
then subsequent iterations could have utilized the calculated 
velocities to determine a second approximation for streamline 
curvature. However this approach was considered excessively 
time-consuming for the present problem and not sufficiently 
justified by the requirements of computing vented Coanda flows 
where the nozzle spacing is not large and the wall curvature is 
substantial. The major curvature effects are experienced in a 
region close to the wall itself. 

4. The Effect of Streamline Curvature on Jet Mixing Turbulent Shear 
Stresses - As pointed out earlier, the first-order effects of 
curvature on turbulent mixing length have been estimated by 
Bradshaw. An approximate expression for the dependence of mix- 
ing length upon Richardson number is included in the calcula- 
tion method. However this approximation is not well validated 
by experimental data and includes curvature effects significant- 
ly larger than those considered by Bradshaw, hence, this adds 
another element of uncertainty to the flow calculation. 

Figure 21 shows the effect of varying the Richardson number 
coefficient (from 3 to 10) upon a computed static pressure distribution along 
the wall for case 6. This range of Richardson coefficient may be thought to 
represent the uncertainty in the magnitude of the effect but suggests very 
little alteration on computed wall static pressure distributions. Increasing 
the value of the Richardson number coefficient tends to decrease the turbulent 
shear stresses in the upper part of the jet mixing zone and to increase them 
on the lower side. 

Figure 22 shows the sensitivity of the calculation of wall static 
pressures for case 6 upon the assumed value of the streamline curvature decay 
coefficient used in equation 39. The effects of variation of this coefficient 
are naturally unimportant in the downstream region where curvatures are small 
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but can be quite significant in the region for a less than 0, i.e., close to 
the zones of high wall curvature. 

5 . 3 Loc us of Maximum Stagnation Pressure 

Figures 23 through 30 show the computed location of the line of 
maximum stagnation pressure from the nozzle to a point far downstream in the 
test section. Also shown are test data points taken from the maximum stagna- 
tion pressure in the upstream zone and from the maximum velocity, i.e., maxi- 
mum stagnation pressure for the downstream region in which curvature effects 
are negligible. 

As with the static pressure comparisons the best agreement between 
calculated and measured values of the locations of maximum stagnation pressure 
correspond to those experimental cases in which there was the largest degree 
of Coanda turning and the smallest spacing between the nozzle and the wall. 

This is shown particularly by the comparison for cases 7 and 8. In other cases, 
for example, cases 5 and 6 (with large venting between the large spacing be- 
tween the nozzle and the wall) the computed location of maximum stagnation 
pressure shows a substantial deviation from the experimental results. These 
two cases as pointed out earlier appear to show a substantially diminished 
Coanda effect. The jet clearly does not cling as closely to the wall as the 
computer model predicts. In general the differences between computed and ex- 
perimental results may be ascribed to the reasons mentioned earlier for the 
static pressure discrepancies. 

, ,, Figures 31 and 32 show the total pressure profiles across the mix- 

ing section. The analytical prediction shown with continuous lines shows a 
small deviation from the experimental results in the tertiary flow that is 
near the bottom wall. 

5 . 4 Velocity Profiles 

Figures 33 through 40 show comparisons of non-dimensional velocity 
profiles at various locations throughout the test section for each of the eight 
cases investigated. Owing to the uncertainty in velocity determination as 


evidenced by the Integrated mass flow of discrepancy being up to 15% the large 
uncertainty level must be attached to each velocity determination. Hence, the 
rather substantial discrepancies between experimental and computed velocities 
are not conclusive indications of the degree of reliability of the analytical 
method. The experimental velocity profiles for Runs 5 and 6 show that the maxi- 
mum velocity region has been shifted towards the bottom wall indicating a re- 
duced coanda effect at large nozzle spacing. 


36 



Section 6 


GENERAL CONCLUSIONS 

1. An approximate method has been developed for calculation of vented Coanda 
flows in ducts. A method has been confirmed by experimental data with an- 
gles of turning up to 67.5° and for close spacing between the nozzle and 
the curved wall. At larger spacing the model indicated flow separation 
which limits the availability of the model to represent the flow profile 
in the downstream zone. 

2. A quasi-orthogonal method of computation, which is more rapid than an iter- 
ative solution of the elliptic boundary value problem, appears best suited 
to ducted Coanda flows with low venting and large curvature. It requires 
approximate specification of streamline curvature decay with distance from 
the wall, and thus is best suited to cases in which the jet sheet is lo- 
cated close to the wall. It is desirable to extend the use of the method 
to non-vented Coanda flows. 

3. Though the effects of streamline curvature on mixing length are known only 
for small curvature, and perhaps uncertain within a factor of 3, a simple 
correction for mixing length in terms of Richardson number appears to pro- 
vide a reasonable estimate for jet curvatures d/R of the order of 0.02. 

4. The flow model developed provides good agreement with secondary and pri- 
mary mass flows measured with orifice plates. Integration of velocity pro- 
files failed to provide satisfactory agreement with orifice measurements of 
mass flow apparently due to the effects of a high turbulence and high shear 
in the mixing zone on the stagnation pressure readings of a Kiel probe, 

5. The flow model predictions were in good agreement with measured wall static 
pressures except in the immediate region of the nozzle apparently due to 
upstream boundary layer and nozzle thickness effects, and due to incipient 
flow separation in certain of the tests. 
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A sensitive measure of the degree of agreement between flow model and ex- 
perimental results is the location of the maximum stagnation pressure line 
in these highly curved flows. 



APPENDIX A 

Sample Tabulation of Static Pressures (Tests 7 and 8) 


Tap Position 
inch 

- 5.50 “ 

- 5.00 

- 4.50 

- 4.00 

- 3.50 

- 3.00 w 

- 2.50 

- 2.00 

- 1.50 

- 1.00 

- 0.50 

0.00 ^ 
0.50 
1.00 

1.50 

2.00 

2.50 

3.00 

3.50 

4.00 

4.50 

5.00 

5.50 

6.50 

7.50 

8.50 

10.00 

11.50 

12.50 

14.50 


Wall Static Pressure 
Run 7 


of water gage 


Top Wall 


Bottom Wall 


- 8.82 
-10.25 
-9.80 

- 5.75 

- 2.85 
+ 2.85 
+ 7.00 
+ 9.75 
+11.90 
+T3.65 


9.35 

9.55 

6.00 

2.55 
+ 2.85 
+ 7.00 
+ 9.75 
+11.90 
+ 13.65 


Rui 

1 8 

Top Wall 

Bottom Wall 

- 

- 2.10 

- 

- 2.60 

-12.30 

- 2.95 

-13.70 

- 3.65 

-16.10 

- 4.15 

-18.95 

- 4.30 

-21.30 

- 4.63 

-22.95 

- 4.90 

-24.30 

- 5.85 

-25.50 

- 6.95 

-28.97 

- 8.45 

-36.99 

-10 . 45 

-33.76 

-13.15 


-15.80 

-19.80 

-32.30 

-23.30 

-31.96 

-25.70 

-30.06 

-26.59 

- 

-28.29 

-32.44 

-29.17 

- 

-30.23 

-32.30 

-30.94 

-32.16 

-31.62 

-33.15 

-32.91 

-37.06 

-35.43 

-37.06 

-36.04 

- 

-36.96 

-31.48 

-31.42 

-26.59 

-25.88 

-17.85 

-17.60 

- 

-11.18 

- 7.65 

- 7.05 

- 4.10 

- 3.75 

- 1.95 

- 1.30 












APPENDIX B 

Finite Difference Equations 

This Appendix provides the detailed derivations of the finite dif- 
ference equivalents of the momentum and energy conservation equations (5) and 
(7) respectively. For convenience the following definitions are introduced: 

S = pu(y + pE) 

and 


S’ 




These definitions permit the momentum and energy equations to be expressed 
as : 


3u _ _ dP , _3^ 

3X “ 2p dX “ 3i|; [f 3i^ 


IT ^ Y - 1 ^ Y - 1 

3X 2pY ^ dX Y 


uS 



+ u 



(B-1) 

(B-2) 


Before approximating these equations with finite difference rela- 
tions a system of grid lines parallel to the X and ip axes must be introduced. 
As illustrated in Figure B-1, a nodal point coincides with each intersection 
of these lines. Lines parallel to the ip axis are termed m-lines and those 
parallel to X axis n-lines. Each node is given a double subscript, the first 
being the number of the m-line passing through it, and the second the n-line 
number. 

m=l m=2 m=»3 



Figure B- 1 Definition of Grid Lines for Finite Difference Solutton 
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The values of the variables on the m = 1 line are the known initial conditions. 
The conservation equations express for each node on the m = 2 line its inter- 
relation with other nodes on the m = 2 line and nodes on the m = 1 line. If 
m = 2 line nodes are only related to nodes which lie on the m = 1 line, the 
finite difference scheme is termed explicit. If an m = 2 node is also related 
to a number of other m = 2 nodes, the scheme is termed implicit (See Figure 
B-2) . 

? ^ 

\ 

} 

t 

> -J> 

IKPLICIT 
Figure B-2 

Diagrams of Explicit and Implicit Solutions 

The implicit form of finite difference schemes leads to a series of N simul- 
taneous algebraic equations relating the known initial conditions on the m == 1 
line and the unknown variables on each of the N nodes on the m = 2 line. After 
solution of these simultaneous equations, the variables on the m = 3 line are 
expressed in terms of the known values on the m = 2 line. Proceeding in this 
manner, a solution to the complete flow field is marched out. Although simpler 
to program, the explicit scheme shows unstable characteristics if the m-lines 
are widely spaced relative to the n-line spacing. Implicit schemes show much 
more stable characteristics and therefore allow much larger m-line spacings, 
thus reducing computation times. The computer procedure presented in this 
report employs a system of implicit finite difference approximations which are 
defined using the notation described in Figure B-3. 
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— m-lines 


n-1 — 


n-lines 


Figure B-3 

Implicit Finite Difference Term Definition 

The velocity at nodes n + 1 and n-1 can be expressed in terms of a Taylor 
Series expanded about node n, on the sa:me m-line. 


n + 1 


9u 

u + — 

n 1 


a^u 

+ + higher order terms 

n ^ n 


)" 2 

u = u - Ai|;_ — + — ^ + higher order terms 

n-1 n 2 d\p 2 _,2 


Combining these equations to eliminate — ^ yields, 




(A^l) 


= -f +1^ I (AD-^A.>2^ + Atlj^Aip^") 


n + 1 


+ higher order terms 




Neglecting terms of the order (A^) and higher, yields 


3u 

8^ 


“n + 1 [A 4 I 2 J - 1 [aiIIj^ Aifi^J "n 


n 


Ai/^2 + 




®5 ■ A<.2(A'('2 + AV 


and S 


All/, 


4 Ai|>j^(A\|»2 + 


yields 


iM 

di|; 


®4 (“n + 1 ■ “n^ ®5 <% ‘ "n - 1> 


(B-5) 


Similarly, 


If -S4<^n + l- 

^ n 


(B-6) 


The second derivative term in the momentum equation is approximated using the 
following Taylor Series expansions, 


fsiu] „ sM 


f3ui . J. 

84> 


S — 

I 


n 




4 3t' 


S — 

r sf'j 


+ higher order terms 
Ai|; 


(B-7) 


fs i“l * fs 1 h1 2 -L [s 


^*1 3^ 


4 A .2 


s — 1 

I 


n 


+ higher order terms 


(B-8) 
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Neglecting terms of the order of and higher yields. 


_L 

3if; 





s — 

n 



S —1 

+ 1/2 


9u 

3i|; 


n - 1/2 






(S +S)(u .n-u) 

n + 1 n n + 1 n 


Alp. 


(S + S .)(u - u .) 

n n-1 n n-1 

Aip^ 


(B-9) 


Similarly, 


9 

c. 9T 

1 

(S' + S' )(T ^ - T )1 

n + 1 n n + 1 n 

9ip 

f 9t 

L ' *'^1 + ^'^2 

Aip^ 


(S' + S’ 

n n - 



T 

n 



(B-10) 


The velocity at a node located at the intersection of the down- 
stream m-line and any n-line ^ can be expressed in terms of the following 

Taylor Series, 


u 


2, n 


^1, n 9X 


AX + 


9^u 


n 


9X 


(AX) + higher order terms (B-11) 


Use of the boundary layer equations implies that gradients in the X- direction 
are much smaller than those in the ip- direction. Therefore it is permissible 
to use a simpler approximation of the X- direction derivatives. 

2 

Neglecting terms of (AX) and higher yields. 


9u | _ ^2, n ^1, n 

9X AX 

'n 

This approximation is termed "backward-difference". 


(B-12) 
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similarly 


I T - T 
^ 2« n 1, n 

9X AX 

'n 


(B-13) 


The only terms in the energy and momentum equations which cannot be approximated 
using the preceding equations are those containing the pressure gradient dP/dX. 
Assuming this gradient varies linearly throughout the AX interval yields: 


dX 


dP 


dX 


m 


+ 


dX 


m 


(B-IA) 


Momentum Equation 


Combining equations (B-1) , (B-9) , (B-12) and (B-14) yields 


u 


1, n 


^"2. n - "l. 
AX 


4p 


1, n 


dX 


m - 1 


dP 

dX 


"I 


m = 2 


+ 


‘1 , n 
2^. 




(S^ . . + S^)(u_ „ . 1 - u. (S„ + ^)(u_ ^ - u. ^ 

n + 1 n 2, n + 1 2 , n _ n n - 1 2 , n 2, n - 1 


Alp. 


Alp, 


(B-15) 


This equation can be expressed in the form 


A ,u + B u^ 1 + C -u_ 1 = B 1 

n~12, n n-12, n + 1 n~12, n~l n-1 


(B-16) 


in which the coefficients are defined by equations (23) through (28) of the main 
text. 


Energy Equation 

Combining equations (B-2) , (B-5) , (B-10) , (B-13) and (B-14) yields 
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« 

1 > n Zy n 


AX 


T ) 

1, nV. Y - 1 

= Ut S_ 

Y 1, n 1, n 


« 4 - 1 ~ ) 

4 z, n + 1 z, n 


n “ “ 2 , n - 


+ u 


1, n 


+ AiJ^2 


(S' . + S’ )(T„ . T - T. ) 

n + 1 n 2, n + 1 2, n 




^’n r, " ^9 r, 

n n — 1 z , n z , n - 1 


n rY - 1 


Alp, 


u. 


JC 1> n 


4p 


1, n 



dP 

4- ^ 


dX 

T ^ dX 

J 

_ 

m = 1 


m = 2 


(B-17) 


This equation can be expressed in the form: 


n - 1 


2, n 


+ B 


n - 1 


2, n + 1 


+ C 


n - 1 


2, n - 1 


= D 


n 


- 1 


(B-18) 


in which the coefficients are defined by equations (31) through (36) of the 
main text. 
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APPENDIX C 
Solution Procedure 

The calculation procedure starts at the upstream flow boundary, 
where the values of all flow variables must be known or assumed. Specifica- 
tion of the velocity and temperature distribution, dimensionless eddy viscosity, 
duct and nozzle inlet dimensions, and working fluid, defines all initial con- 
ditions . 


The known initial conditions, m = 1 line, are related to the un- 
known conditions, m = 2 line, by the previously derived equations, and assumed 
boundary conditions. These inter-relations form a set of n - 2 simultaneous 
algebraic equations, where n is the number of n-lines, and the equations are 
shown in Appendix B. The resultant matrix of coefficients is tridiagonal in 
form except for the initial and final rows which only contain two terms. Rapid, 
exact solutions to this type of matrix are obtained using the Thomas Algorithm, 
a successive elimination technique, which is described in this Appendix. 

The solution for the variables on the m = 2 line is iterative, be- 
cause of the presence of the unknown pressure in the momentum equation. The 
procedure adopted was to estimate the pressure gradient, and solve the equations, 
using the algorithm. The equations automatically satisfy conservation of mass, 
momentum, and energy, but only one pressure gradient yields the correct wall 
geometry. The duct dimension corresponding to the estimated pressure gradient 
was calculated from the m = 2 line variables. The pressure gradient was then 
incremented by a small percentage of its initial estimated value, and the cal- 
culation process repeated for a new duct dimension. A third estimate of the 
pressure gradient was obtained by interpolation between the two calculated, 
and the actual duct dimension. In almost all the calculations performed to 
date^ this value has been acceptably close, within 0.001%, to the actual duct 
dimension. If this criterion is not met, a further iteration is applied, and 
a fourth solution obtained. 
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The now known variables on the m = 2 line become the new m = 1 line 
variables and the procedure is repeated for another set of m = 2 line variables. 
Thus a solution to the complete flow field is marched out. 


The difference form of the momentum and energy equation is: 


A .X + B .X . , + C = D , 

n-ln n-ln + 1 n-ln-1 n-1 


(C-1) 


where X is either u or T. If the number of n-lines is n, there are n - 2 
equations of the form (1) and two equations expressing the boundary conditions. 
The first and the last equations represent the boundary conditions, which are: 


and 


u, = 


u = 0 
n 





= 0 


n 


(C-2) 


(C-3) 


Equation (C-2) can be written in terms of X as follows: 

X^ = X^ = 0 (C-4) 

Equation (C-3) correspondingly becomes: 

= X^ (C-5) 

and 

X = X (C-6) 

n-1 n 

Equations (C-4), (C-5) and (C-6) can be written in terms of X as follows: 


X = KX , (C-7) 

n n-1 

X^ = KX^ (C-8) 

where K is 0 for the momentum equation and unity for the energy equation. 

Thus, the matrix form of the equation (C-1) is shown on the following page 
(Table C-1). 
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Table C-1 


Matrix Form of Equation C-1 Designated as Equation C-1 


1 — 











"~r 

r n 



'1 

-K 

0 

0 

0 

- 

0 

0 

0 

0 

0 

0 



0 




0 

0 

- 

0 

0 

0 

0 

0 

0 

^2 



0 


^2 


0 

- 

0 

0 

0 

0 

0 

0 



°2 

0 

0 

s 



- 

0 

0 

0 

0 

0 

0 



°3 

- 

- 

- 

- 

- 

- 

- 

- 

- 

- 

- 

- 

Vl 

= 

V2 

0 

0 

0 

0 

0 

- 

c 1 

n-1 

^ 1 
n-1 

Vi- 

0 

0 

0 

X 

n 


Vl 

- 

- 

- 

- 

- 

- 

- 

- 

- 

- 

- 

- 

X .1 

n+1 


D 

n 

0 

0 

0 

0 

0 

- 

0 

0 

0 

C o 
n-2 

A 9 

n-2 

S 9 

n-2 

X 1 

n-1 


V2 

0 

0 

0 

0 

0 

- 

0 

0 

0 

0 

-K 

1 

X 

__n _ 


0 



The second equation Is: 


CiXi + + 3^X3 = 


(C- 9 ) 


Substituting equation (C- 4 ) into this equation yields: 


A- 3X2 + 3^X3 = 


(C- 10 ) 


where 

th 

The n*^ - 1 equation is: 


C .X _+A _X ,+B „X=D - 
n-2n-2 n-2n-l n-2n n-2 


(C- 11 ) 


Substituting equation (C- 7 ) into this equation yields; 


C .X + X = D 

n-2n-2 n-2n-l n*~2 


(Crl 2 ) 


where A'=A»+KB 

n-2 n-2 n-2 

Thus the n equations (C-8) can be reduced to the n -2 equations shown on 
Table C- 2 . 

The Thomas Algorithm 


Starting with the first equation, X^ can be expressed in terms of 
X„. The second equation gives X- in terms of X,. Continuing through all the 
equations until the n -3 equation gives X^ 2 terms of 2* Com- 

bining this with the last equation gives _ 2* Working backwards through 
the equations then allows the remaining unknowns to be found. This procedure 
is most easily applied by defining the following: 

D, 


”1 = ^'1 


81 = 



- 1 

n - 1 
n - 1 


n = 

2, 

3 - - - (n - 

2) 

(C- 14 ) 

w 

n 

= A 

n 

- C Q 
n n - 

1 

n - 

2, 

3 (n - 

2) 


®n 

= D 

n 

c g 

n n - 
W 


n = 

2, 

3 (n - 

2) 



n 
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Table C-2 


! 

1 


Matrix Form of Equation C-8 with Simplified Terms Designated as Equation C-13 




0 

0 

- 

0 

0 

0 

0 

0 

0 


r n 






0 


0 

0 

0 

0 

0 

0 


^3 


■>2 

0 

"3 

^3 

®3 

- 

0 

0 

0 

0 

0 

0 





- 
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- 

- 

- 

- 

- 

- 

- 

- 

- 


Vi 

= 

“n-2 

0 

0 

0 

0 

- 

C 1 

n-1 

^ 1 
n-1 

B _ - 
n-1 

0 

0 

0 


X 

n 


D T 
n-1 



- 


- 

- 

- 




- 


\+l 


D 

n 

0 

0 

0 

J 0 

•- 

0 

0 

0 

n-3 

n-3 

n-3 


\-2 


D o 
n-3 

0 

0 

0 

: 0 


0 

0 

0 

0 

C 9 

n-2 

n-2_ 


Vi_ 


^ 9 

_n-^ 




I 


Ul 

!-» 



Equations (C-13) then reduce to: 

X = g and X = g . - Q„ -X n = (n - 2) , (n - 3) , 2 

n-1 n~2 n n-1 n-ln + 1 

(C-15) 

If the values of W, Q and g are calculated in order of increasing n using 
euqations (C~14) , then equations (C-15) can be used to calculate the values of 
X in order of decreasing X starting with \ To clarify this procedure, 

the method is now used to solve the following four simultaneous equations; 


S ""2 ®2 ° 

° *^3 S ®3 

0 0 A-^ 



^2 






°2 


""a 


^3 


3_ 


”4 


A'i^2 + 


Wi = A>^ 

Q 

^1 


hence ^2 ^1 ~ ^1^3 


V3 ^ ®2^4 ^ "2^2 = ^2 


'^2 = "^2 - '^ 2'^1 
^2 

"^2 - W 2 


S 2 


^2 ~ *^2^1 


hence X^ = g^ - X^p2 


(C-16) 
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^ ^ S''3 = ®3 


^2 *3 ” *^ 3*^2 


<>3 4 




^3 “ ^ 3^2 

W. 


hence X^ = g^ - Q^X^ 


A',X- + C,X, = D, 
4 5 4 4 4 


^ = ^4 - "=4S 
“4 - '^483 


(C-17) 


84 ■ 

hence X - g, (C-18) 

5 4 

Substituting in equation (C— 16) yields X^. Equations (C— 17) and (C— 18) are 
special forms of equations (C-15) for n = 6 and n = 4, 
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Appendix D 


COMPUTER PROGRAM 


The computer program is designed to analyze the flow in a duct 
where geometry is shoxra in Figure D-1. 



The computation proceeds from the initial line by moving along the 
lower and upper wall a specified distance and then defining two arcs from the 
new wall locations to the duct midpoint. Since only a rough approximation Is 
available for the wall slopes the two arcs do not necessarily meet resulting 
in the wall points becoming slightly unsynchronized as the computation pro- 
ceeds . 
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The program organization consists of a main program (NAS) and 
sixteen subroutines and function subroutines. The functions of the main pro- 
gram and subroutines are: 

PROGRAM NAS 


The main program is divided as follows: 

Input and Data Initialization: Cards $NAlOto $NA2200 

This section initializes the constants of the program, reads and 
prints the computation conditions and duct geometry, and puts this data in 
nondimensional form. 

Initial Conditions; Cards $NA2210 to $NA2810 


The subroutine INCOND is called to define the starting conditions. 
The initial flow conditions are then put in dimensional form and printed. 

Main Body of Program: Cards $NA2820 to $NA4810 

The computation proceeds down the duct in a sequence of steps. 
Values of pressure, temperature, velocity, density etc. are computed which are 
consistent with the previous step values and the geometry of the duct. The 
process stops when the end of the duct is reached. 

Eddy Viscosity; Card $NA3Q1Q 

The dimensionless eddy viscosity is calculated in subroutine EDDY 
using data from the preceding step. 

Streamline Step Size: Cards $NA3100 to $NA3400 

An appropriate step size along the duct is determined for the 
streamline of maximum velocity. Consistent step sizes are then 
determined for all other streamlines. 
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Pressure Gradient Approximation: Cards $NA3730 to $NA4540 


The pressure gradient is determined by selecting a value such that 
the computed duct widths minus the actual duct width equals 0 + 
. 0001 . 


Flow Boundaries: Cards $NA4550 to $NA4700 

At each step the boundaries of the flow regions are checked to 
determine when shear layers vanish. 

Output Section; Cards $NA4820 to $NA5750 

The flow variables are presented in dimensional form at preselected 

intervals . 


Figure D-2 is a flow chart of the main progran NAS. 
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SUBROUTINES 


INCOND 

This subroutine contains the computation which defines the flow 
conditions at the initial station (Fig. D-3) . The steps used in the process 
are: 


1. Points on the boundaries (x^, y^) , (x^, y^) are defined such 

that two circular arcs are normal to the wall at (x^^, y^) and 
(x^, y 2 > and pass through the nozzle ^noz^ * 

2. A flow split is determined. 

3. The subroutine OMGSET is called to set the streamline loca- 
tions. 

4. The subroutine TMPSET is called to determine the temperature 
distribution. 

5. The remaining flow variables are calculated. 

6. The location of the nozzle streamline checked to see if it is 
within tolerance of its specified location. If it is then the 
computation is returned, if not a new flow split is determined 
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OMGSET 


This subroutine sets the distribution of streamlines. The flow is 
divided into regions of width DOMGl, D0MG2, D0MG3 (Fig. D-4 and D-5) . Given 
an initial spacing at each flow edge (DMHALL and DMJET) and a specific number 
of points (KWALL, KJET) the subroutine fills in the spaces close to the edges. 


OMGN 

DMWALL = KWALL 


0MGJ2 DHJET = KJET 

DMJET = KJET 


0MGJ3 J>MJET P KJET 

DMJET E KJET 


0.0 


DMWALL = KWALL 


Figure D-4 
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TMPSET 


selected. 



Figure D-5 
Flow Structure 


An Initial value of static temperature at the nozzle, TS(NMID), is 
Values of the temperature are then computed at the outer wall con- 
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sis tent with the wall curvature. The distance between the wall and nozzle 
consistent with TS(NMID) is returned to the calling program. 

WALLS 


Wall curvatures are calculated from the wall geometry data. This 
data is then used to interpolate values of wall location slope and curvature 
for specific axial locations. The computer calculates curvatures assuming 
the data is segmented with at least three points per segment. When five or 
more points are available a least squares parabola using the data points and 
the two points to either side is used to smooth the data. Subsequent inter- 
polations are parabolic. 

LSQ 


This subroutine calculates the least squares parabolas for walls. 


EDDY 


This subroutine contains the eddy viscosity calculation. 


DSMOVE 


This subroutine contains the calculation of the distance moved 
along the walls in each step. 

SDWE 


This subroutine contains the calculation of temperatures and 
velocities at a computing station. 

CALC 


This is a simultaneous linear equation solution subroutine. 


CORDS 


The distance between streamlines and the x, y coordinates of the 
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streamlines are computed in this subroutine. The value DDIST, the difference 
between the y coordinates of the two arcs at the middle duct point, is returned, 

LOOK 


The algorithm identifies an edge of a core region and then checks 
to see if a core region that existed in the previous step has disappeared. If 
this has occurred that core width is set to 0.0. Core regions may disappear 
in any order. 

REMOVE 


This subroutine is used to remove grid points as the computation 
proceeds in order to reduce computing time. This is accomplished by every 
tenth step scanning the velocity array to see if the velocity gradients in 
the shear regions are less than a tolerance. 

ARCDIS 


This subroutine defines wall coordinates such that the arc is nor- 
mal to the wall and nozzle. 

FUNCTION SUBROUTINES 

CENTRE 


This is used in conjunction with ARCDIS. Given a nozzle and wall 
Ideation it projects the two tangents, finds their intersection and returns 
the difference in length. 

CURVDS 


Used in DSMOVE to calculate distances along an arc. 


SERIES 


Used in OMGSET to determine streamline spacing in the boundary 
and shear layers. 
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INPUT DATA 


Card 1 - Format I1,18A4 

NUNIT ~ units Indicator 

NUNIT = 0 English units 
NUNIT ^ 0 S.I. units 

TITLE(I) - title, printed at start of first output page 
Card 2 - Format 6F10.0 


POl - primary stagnation pressure - psia or pascal 
TOl - primary stagnation temperature - Rankine or Kelvin 
P02 - secondary stagnation pressure - psia or pascal 
T02 - secondary stagnation temperature - Rankine or Kelvin 
MASSl - primary mass flow rate - Ibm/sec-in or kg/sec-m 
MASS2 - secondary + tertiary flow - Ibm/sec-in or kg/sec-m 

Card 3 - Format 515 

K1 - no, of grid points in primary flow (even) 

K2 - no. of grid points in secondary flow (even) 

K3 - no. of grid points in tertiary flow (even) 

NPCYCL - print cycle, (e.g. NPCYCL = 10 causes print every ten 
steps) 

NQUICK - 0 - full print; 1 = partial printout 
Card 4 - Format 4F1Q.0 

XNOZ - X coordinate of nozzle center - in. or meters 
YNOZ - y coordinate of nozzle center - in. or meters 
TNOZ - nozzle angle - degrees 

N - nozzle slot width - in. or meters 

Card 5 - Format 1115 

MST - number of segments used to describe lower wall geometry. 


63 


10 2: mst 

NDI(I, 1) - MST values of the data point number of the segment end 
point, NDI(MST,1) = NPAIRl, the total number of points 
needed to describe the wall. 

NDI(I + 1,1) - NDId, 1) ^ 2 

Cards 6 - Format 2F10.0 - NPAIRl Cards 

XW(I, 1) - X lower wall coordinate 
YW(I, 1) - y lower wall coordinate 

Card 7 - Format 1115 

MSS - number of segments used to describe upper wall geometry 
10 ^ MSS ^ 1 

NDI(I,2) - MSS values of the data point number of the segment end 
point, NDI(MSS, 1) = NPAIR2 the total number of points 
needed to describe the wall. 

NDI(I + 1,2) - NDI(I, 2) 2 

Cards 8 - Format 2F10.0 - NPA1R2 Cards 


XW(I, 2) - X upper wall coordinate 
YW(I, 2) - y upper wall coordinate 

Card 9 - Format I5,1QF5.0 


NFULL - number of values of x for which full output is required 

XPR0F(I) - NFULL values of x for which full output is required 

The input data Kl, K2, and K3 (Card 3) denote the number of grid 
points in the primary, secondary and tertiary flow. These are arbitrary num- 
bers chosen so as to give the desired spacing of output data in the three re- 
gions. A total Kl + K2 + K3 of between 20 and 30 should give good results. 

The data numbers MST, MSS (Cards 5 and 7) denote the number of 
segments needed to describe the wall geometry. For walls with continuous 
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curvature values, a value of 1 Is sufficient. Values of MST and MSS greater 
than 1 allow the user to describe walls with discontinuous slopes and curvatures. 
The data NDI(1,2) and NDI(I,1) are the data point numbers at the boundary seg- 
ment ends. 


A san^le of Input data Is shown In Table D-1. 
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Table D-1 


1 

TEST CASE 

253934. 

322. 

0 

10 10 

-.07366 

.04023 

6 

12 16 

-.1524 

-.10396 

-.1511 

-.1024 

-.1499 

-.10033 

-.1397 

-.00694 

-.1270 

-.07450 

-.1143 

-.06535 

-.1016 

-.05837 

-.0889 

-.05339 

-.0762 

-.05004 

-.0635 

-.04712 

-.0597 

-.04678 

-.0587 

-.04671 

-.0500 

-.04517 

-.0254 

-.04023 

0.0000 

-.03523 

.01905 

-.03160 

.02159 

-.03115 

.02032 

-.03137 

.02549 

-.03048 

.03810 

-.02870 

.05000 

-.02736 

.06350 

-.02647 

.07620 

-.02596 

.08636 

-.02564 

.08763 

-.02560 

.08890 

-.02558 

.10160 

-.02536 

.12700 

-.02492 

.15240 

-.02448 

.17780 

-.02404 

.19050 

-.02383 

.20320 

-.02383 

.22860 

-.02383 

.25400 

-.02383 

.26670 

-.02383 

.27940 

-.02449 

.30480 

-.02582 

.35560 

-.02848 

.40640 

-.03115 

.45720 

-.03381 

.50800 

-.03649 

.55080 

-.03913 

.58420 

-.04038 


Sample of Input Data 
5 METRIC 

301. 

5 0 

-35. .0030ft61 

36 31 35 43 


1.696 


7.393 
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Table D-1 


5 

16 26 

«.0845 

.156410 

-.0838 

.143230 

-.0826 

.134200 

-.0794 

.120680 

-.0762 

.111240 

-.0699 

.097190 

-.0635 

.086450 

-.0572 

.077670 

-.0508 

.070240 

-.0381 

.058260 

-.0254 

.049090 

-.0127 

.0^2030 

o.oooo 

.036730 

.01270 

.032960 

.01588 

.032240 

.01905 

.03160 

.02032 

.03137 

.02159 

.03115 

.02549 

.03048 

.03810 

.02870 

.05080 

.02736 

.06350 

.02647 

.07620 

.02596 

•08636 

.02564 

.08763 

.02560 

.08890 

.02558 

.10160 

.02536 

.12700 

.02492 

.15240 

.02448 

.17780 

.02404 

.19050 

.02383 

.20320 

.02383 

.22860 

.02383 

.25400 

.02383 

.26670 

.02383 

.27940 

.02449 

.30480 

.02582 

.35560 

•02848 

.40640 

.03115 

.45720 

.03381 

.50800 

.03649 

.55880 

.03913 

.58420 

•04038 


1 .038 


31 


(Concluded) 
35 43 



OUTPUT DATA 


The program printed output consists of: 

1. Title 

2. Input flow and geometry (Cards 2, 3, and 4) 

3. Wall Geometry 

X, Y - wall coordinates (smoothed) 

CURV - negative of wall curvature 

4. Values of X at which full output is required (values of 
XPROF(I) from Card 9) 

5. Flow split 

SPLIT “ ratio of tertiary flow to the sum of secondary and 
tertiary flow 

6. Initial conditions along the computing station through the 
nozzle 

ISENTROPIC NOZZLE THRUST PER UNIT WIDTH - lb /in or N/m 

AMASSl, AMA.SS2, AMSS3 -* primary, secondary and tertiary flow rates - Ib/s/in 

or kg/s/m 
J “ station number 

X, Y - coordinates of computing station streamline inter- 
section - in or m 

DN - distance from lower wall - in or m 
THETA - streamline angle - degrees 

K - negative of streamline curvature - 1/in or 1/m 
OMG - streamline function 

PO, PS - total and static pressure - psi or pascal 

TS - static temperature - degrees Rankine or Kelvin 
RHO - density - lb/£**3 or kg/m**3 
U - speed - f/s or m/s 
M - Mach number 
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7. Flow description downstream - partial output 

NSTEP - step number 

XI (XO) - coordinate along lower (upper) wall - in or m 

PI(PO) - static pressure at lower (upper) wall - in of 
H 2 O or pascal 

USTARI(USTAR0) - friction velocity at lower (upper) wall “ f/s 
or m/s 

KI(KO) - negative of wall curvature at XI (XO) - 1/in or 
1/m 

RNI(RN0) - Richardson number at lower (upper) wall 
UMAX - maximum velocity, f/s or m/s 
NUMBER OF ITERATIONS 

DELTA X - step size at lower wall - in or m 
DELTAS I ( DELTAS 0) - distance increment of lower (upper) wall - in 

or m 

SI (SO) - cumulative distance along lower (upper) wall - 
in or m 

SHRIN(SHROR) - shear stress at lower (upper) wall - psi or 
pascal 

TOTAL AXIAL MOMENTUM PER UNIT WIDTH - Ibf/in or N/m 

THRUST RATIO - 

8. Flow description downstream - full output 
Same as partial output plus 

N1 . . . N6 - streamline numbers at edges of boundary and 
shear layers 

COREl, C0RE2, C0RE3 - widths of primary, secondary, and tertiary 

regions - in or m 

DELI, DEL2, DEL3, DEL4 - widths of boundary and shear layers - in or m 

J - streamline number 
1/R - negative of curvature - 1/in or 1/m 
THETA - angle of streamline - degrees 
X, Y - computing mode location - in or m 
YREL - relative y coordinate 
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UREL - relative velocity (UMAX is normalizing quan- 
tity) 

POREL - relative total pressure 

PS - static pressure - psia or pascal 
E - Eddy viscosity 

TOTEMP - total temperature - degrees Rankine or Kelvin 

In addition several warnings may be printed. 

NO CONVERGENCE, DPDSA = ... DPDSB = ... DDISTB ... STEP 4 

Convergence was not achieved in establishing the pressure gra- 
dient. The criteria for convergence is that [ DDISTB | the dis- 
tance between the computed and actual will location be <.0001. 
DPDSA and DPDSB are pressure gradient increments. 

NEGATIVE SHEAR STRESS AT THIS STATION. THIS INDICATES POS- 
SIBLE SEPARATION. SUBSEQUENT RESULTS SHOULD BE USED WITH CAU- 
TION. 

The message is self-explanatory. When this occurs the calcu- 
lation may be unstable and may stop for a variety of reasons. 
When this occurs the full output at that station is printed. 



•‘^oooooooonooooooooonoonoononooooooonooooonooonoooooooooooooooo 


Figure D-6 


COMPUTER PROGRAM LISTING 


PROGRAM NAS (INPUT • OUTPUT * TAPES* INPUT * T APE6*0UTPUT ) 

THIS PROGRAM CALC»tLATES VEnTED COANOA JET IN A CURVED DUCT 
BY A finite-difference METHOD 

THE PROGRAM CALCULATES THE SPLIT BETWEEN SECONDARY 
AND tertiary mass FLOW 

READ variables 

NUNIT... units INDICATOR— NUNIT*0. -ENGLISH UNITS 

NUNIT GT 0 ..SI UNITS 

poi .. •primary stagnation pressure — psia or pascal 
T01**.PRIMARY stagnation temperature — KANKINE or KELVIN 
P02.. .SECONDARY STAGNATION PRESSURE — PSiA OR PASCAL 
T02.. .SECONDARY STAGNATION TEMPERATURE — RANKINE OR KELVIN 
MASSl ... primary mass FLOW RATE— LBM/SEC-IN, OR KG/SEC-M 
MA5S2... TOTAL OF SECONDARY PLUS TERTIARY FLOW RATE — LBM/SeC-IN. 

OR KG/SEC-M 

K1...NO. OF GRID PTS. In PRIMARY FLOW (EVEN NUMBER) 

K2...NO. OF GRID PTS. IN SECONDARY FLOW (EVEN NUMBER) 

K3-..NO. OF GRID PTS. IN TERTIARY FLOW (EVEN NUMBER) 

NPCYCL... print cycle 

NQUICK... 0 - FULL PRINT* 1 - PARTIAL PRINT OUT 

XNOZ... X COORD OF CENTRE OF NOZZLE ...IN. OR METERS 

YNOZ... Y COORD OF CENTRE OF NOZZLE .. IN. oR METERS 

TNOZ... NOZZLE ANGLE DEGREES 

0... NOZZLE SLOT WIDTH INCHES OR MEtERS 

NPAIRI... NUMPER OF DUCT WALL COORDS* LOWER WALL 
NPAIR2... NUMBER OF DUCT WALL COORDS* UPPER WALL 
XW... array of DUCT WALL X COORDS .. IN. OR METERS 

YW... array of duct WALL Y COORDS .. IN. OR METERS 

NFULL... NUMBER OF PROFILE PRINTOUT LOCATIONS 

XPROF... array of locations WHERE profiles required IN. OR M. 
IT IS ASSUMED THAT P02=PU3* T02=T03 

data variables 

RG...GAS constant — FT.LMF/LBM.R OR N-M/KG-K 
G.. .SPECIFIC HEAT CONSTANT 

NUREF... REFERENCE KIN. VISCOSITY— FT**2/SEC OR M**2/SEC 
TREF... REFERENCE TEMPERATURE— RANKINE OR KELVIN 
RHOREF... REFERENCE DENSl TY— LBM/FT**3 OR KG/M**3 
'GC... CONSTANT— LBM.FT/LrF.SEC**2 OR KG-M/N-SEC**2 
PR...PRANOTL number 
PRT... TURBULENT PRANOTL NUMBER 
TW...WALL TEMPERATURE — RANKINE OR KELVIN 

DIMENSIONLESS VARIABLES 

OMG(I) STREAMLINE FUNCTIONS 

S1*S2»S3*SA*S5(I) OMG DIFFERENCES 
DS(I) STREAMLINE STEP SIZES 

R(I*2) CURVATURES 

THETA(I*2) ANGLES 

PS (I *2) PRESSURES 

U(I*2) VELOCITIES 

TS(I*2) TEMPERATURES 

RHO(I*?> DENSITIES 

VIS(I) VISCOSITIES 

DN(I) DISTANCE TO INNER WALL 

X(t) X COORDINATE 


SNA 0 
SNA 10 
SNA 20 
SNA 30 
SNA 40 
SNA 50 
SNA 60 
SNA 70 
SNA 80 
SNA 90 
SNA 100 
SNA 110 
SNA 120 
SNA 130 
SNA 140 
SNA 150 
SNA 160 
SNA 170 
SNA 180 
SNA 190 
SNA 200 
SNA 210 
SNA 220 
SNA 230 
SNA 240 
SNA 250 
SNA 260 
SNA 270 
SNA 280 
SNA 290 
SNA 300 
SNA 310 
SNA 320 
SNA 330 
SNA 340 
SNA 350 
SNA 360 
SNA 370 
SNA 380 
SNA 390 
SNA 400 
SNA 410 
SNA 420 
SNA 430 
SNA A40 
SNA 450 
SNA 460 
SNA 470 
SNA 480 
SNA 490 
SNA 500 
SNA 510 
SNA 520 
SNA 530 
SNA 540 
SNA 550 
SNA 560 
SNA 570 
SNA 580 
SNA 590 
SNA 600 
SNA 610 
SNA 620 
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c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 


c 

c 


c 

c 

c 

c 

c 


800 

801 

802 

803 

80^ 


9 


YU) 
Md) 
PO (I) 
E(I) 


Y COWHOINATC OF STCAMUINES 
MACM NUMBERS 

stagnation pressures 

EDDY VISCOSITIES 


PRESSURES NORMALIZED BY...P01 
TEMPERATURES, NORMALIZED BY...T01 
VELOCITIES NORMALIZED BY.-.UREFP 
LENGTHS normalized BY ...0 

DYNAMIC VISCOSlTES NORMALIZED BY. • .RHO01*UREFP*D 
KINEMATIC EDDY VISCOSITY NORMALIZED BY...UREFP*D 
STREAM FUNCTION NORMALIZED BY. , .SORT (RHOO 1*UREFP*0) 

MASS Flow rate normalized 8Y..,RH001*UREFP*0 

R(I> IS THE inverse of THE DIMENSIONLESS RADIUS OF CURVATURE 

real NUREF,KSfLM,MA5Sl*MASS2*M( 190) 

DIMENSION XW(99*2) *YW (99»2> »S1 (190) »S2(190) fSSllQO) * S4(I90)* 

1 55U9U) .£(190) *E1 (190) *05(190) *X (190) *Y{ 190) * VIS <1901 * 

2 P0(190) *H(190) *ThETA (190.2) *R ( 1 90 .2) *U U 90 »2) » 

3 TS(19a*?) .RHO(l9o*2) *DIP ( 190) .NX (6) .XPRoF(lO) 

DIMENSION OMG(190) .DN(190),PS(190.2) 

DOUBLE OMG(190) . DN ( 190 ) , PS ( 190 .2) 

DIMENSION PSI (180) 

DOUBLE PSI (190) 

DIMENSION NO! ( 10.2) .CURV(99.2) .TITLEdS) 
ComMON/KURV/MST»mSS*NDI .CIJRV 
CqMMoN/SOLV/ SI »s2*S3.E.EI .H,SA.SS 


IBNA 

SNA 

INA 

SNA 

SNA 

SNA 

SNA 

SNA 

SNA 

SNA 

SNA 

SNA 

SNA 

SNA 

SNA 

SNA 

SNA 

SNA 

SNA 

SNA 

SNA 

SNA 

SNA 

SNA 

SNA 

SNA 

SNA 


Ci30 

640 

650 

660 

670 

680 

690 

700 

710 

720 

730 

740 

750 

760 

770 

780 

790 

800 

810 

820 

830 

840 

850 

860 

870 

eeo 

890 


COMMON/I NEDY/PsI ,RH0, VIS 

COMMON/ 1 NCO/ U,PS»TS»M »P0. theta .R.X.Y.DN 

Common /wall/ npairi* npair2. xw, yw 
COMMON /CONST/ P01.P02.P03. T0l.T02.T03. XI. Xc* YC 
DATA RG.G»GC,TRFF,RHOREF.TW/53.2d .A.32.2.530.,0.075.530.0/ 
DATA NUPEF*pR.PRT/0.0O0l4.0.7.0.9/ 

DATA PRl.PRTl/1. 0.1.0/ 

Data A1»A2»A3/144.. 198.7, ,03611/ 


INPUT data SECTION 


READ (5,800) NUMIT. (TITLE ( D , 1=1 , 18) 
READ(5,801)P01 ,T01*P02,T02»MASS1,MASS2 
READ (5,802) Kl .K2 »k3,NPCYCL.NQUICk 
READ(S, 80i) XNOZ,YNOZ,TNOZ.D 
READ (5,802) M5T, ( ND I (I , 1 ) , 1=1 ,MST ) 
NPAIR1=N0I (MST, 1 ) 

READ(5,6u 3) (XWd,l) ,YW(I,1) ,I = 1,NPAIR1) 

Read (5,802) mss, (ndi d,2 ) ,i=i,mss> 

NPArR2=NDI (MSS, 2) 

read (5,803) (XWd,2) ,YW(I,2) ,I = 1,NPAIR2) 
READ(5,80A) NFuLL. ( XPROF < J) , J=1 ,NFuLL ) 
f ORMAT ( I 1 , 18A4) 

Format (6F10.0) 

FORMAT (1115) 

Format (2fio.O) 


SNA 900 
SNA 910 
SNA 920 
SNA 930 
SNA 940 
SNA 950 
SNA 960 
SNA 970 
SNA 980 
SNA 990 
SNAIOOO 
SNAlOlO 
SNA1020 
SNA1030 
SNAIOAO 
SNA1050 
SNA1060 
SNA1070 
sNAioeo 
SNA1090 
SNAllOO 
SNAlllO 
SNA1120 
SNA1130 
SNA1140 
SNAl 150 
SNA1160 
SNAl 170 


FoRMAT(I 5,10F5.0) SNAllSO 

XPF?0F (NFULL + 1 ) = XW(NPAIR1,1) SNA1190 

WRITE (6,900) (TITLE (I) ,1 = 1,18) SNA 1200 

IF (NUNIT.GT.O) GO TO 9 SNA1210 

WPITE(6,90I) POi ,T01 .P02.T02,MASS1 ♦MASS2.D SNA1220 

WRITE(6,902) XN0Z.YN0Z,TN0Z SNA1230 

GO TO 16 SNA1240 

WRITE(6,914) POI ,T01 .P02,T02,MASS1,MASS2.D SNA1250 

write (6,915) XN0Z.YN0Z,TN0Z SNA1260 

. RG=28A, SNA1270 

Gr=i.. SNA 1 280 
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16 


902 

903 

90^ 


905 

906 
914 


91S 


601 


602 


KI*K2»K3*NPCyCL 


900 

901 


TREF«295. 

rhore:f»i.201 

TWs295. 

NUREFs 1.3E-5 
Alsl. 

A2S110. 

A3S1.0 
WrITE<6*903) 

WRITE(6*904) 

call WALLS<XWf YW« 1000* *0>«T1 tRl vNPAIRI » 1 > 
call walls (XWt YW f 1000. f0.tr2*R2*NPAlR2t 2) 

NPaIRs NPAIRI 

lE CNPAIRI ,6T. NPAjRZ) NPAIR* NPAIR2 
WRITE (6« 905) (J*XW(J«1) tY'^<J*l) «CURV(J»1) >XW(J«?) *YW(J«2) • 

1 CuRV< J*2) • J=l«NPAlR) 

WRITE (6t906) (XPROF< I) ,1*1 ,NFULL) 

FORMATClHl ,39X,43nH*)/40X*lH*f4lX,lH*/40X,lH*,2X* 

1 40HCOANDA EFFECTS IN A PLANE CURVED DUCT */ 

2 40X,lH*,4lX,lH*/40X,43(lH*)//18A4//) 

FoRMAT(25X,29HPRIMARY STAGNATION PRESSURE *,F6,2,5H PSIA/ 

1 25X,26HPRIMARY STAGNATION TEMP. =*F7.2*l6H DEGREES RANKINE 

2/ 25X,31H5EC0NDARY STAGNATION PRESSURE *»F6.2»5H PSIA/ 

3 25X,28HS£C0NDARY STAGNATION TEMP. *,F7«2,l6H DEGREES RANKI 
4NE// 25X,24HPR1 MaRY MASS FLOW RATE *,F8.5,12H LBM/SEC-IN./ 

5 25X,44HT0TAl OF SECONDARY PLUS TERTIARY FLOW RATE =tF8.5* 

612H LBM/SEC-IN.//25X, 19HN0ZZLE SLOT WIDTH =,F6.3,4H IN./) 

FORMAT (25X»21HN07ZLE X COORDINATE *,F10.4,4H IN./ 

Y coordinate =,F10.4,4H IN./ 

ANGLE =»F10.3»8H DEGREES/) 

OF GRID POINTS IN PRIMARY FLOW *,14/ 

OF GRID POINTS IN SECONDARY FLOW =,14/ 

OF GftID POINTS In TERTIARY FLOW =,I4/ 

CYCLE ,14) 


25X,21HN0ZZLE 
25X, 14HN0ZZLE 
FORMAT (25X,39HNUmBER 
2SX»41hnUMBER 
25X,40HNUmBER 
25X, 12HPRINT 


FORMATdHl ,10X,2SHMIXING SECTION DIMENSIONS//// 

1 17X»10HLOWER WALL,35X, lOHUPPER WALL/ 

2 6X»2HJ,10X*1HX,9X1hY, 1 1 X ,4HCURV, 19X, 1 hX ,9X * IHY, 1 IX ,4HCURV 
3/) 

Format ( I7,3X,2F10.4,F15.5.10X,2F10.A,F15.5> 
format (/7X, 12HPROF ILLS AT , 1 OF 10. 5/) 

FORMAT (25X,29HPRIMARY STAGNATION PRESSURE *,El2.5*7H PASCAL/ 

1 25X,26HPKIMARY STAGNATION TEMP. =,F7.2»15H DEGREES KELVIN/ 

2 25X,31HS£C0N0ARY STAGNATION PRESSURE =,E12’.5,7h PASCAL/ 

3 25X,28H5EC0NDARY STAGNATION TeMP. =,F7.2,15h DEGREES KELYI 

4N// 25X.24HPR1MARY MASS FLOW RATE *,F8.5, 9H kG/SEC-M/ 

5 25X,44HT0TAL OF SECONDARY PLUS TERTIARY FLOW RATE =,F8.5, 

69H KG/SEC-M //?5X»19HN0ZZL£ SLOT WIDTH =,F9.6,7H METERS/) 
F0RMAT(25X,21HN0ZZlE X COORDINATE =*F10.6,7H METERS/ 

1 25X,21HN0ZZLE Y COORDINATE =*F10.6*7H METERS/ 

2 25X,14HN0ZZLE ANGLE =*F10.3,8H DEGREES/) 
xnoz* XNOZ/D 

YNOZ* YNOZ/0 

TnOZ* TnOZ / 180.0 * 3.1416 
DO 601 J=1*NPAIR1 
XW(J,1)= XW(J,1) / 0 
YWU.l )* YWC J,l) / 0 
CURV ( J, 1 )=CURV ( J, 1 >*D 
Continue 

Do 602 J*1*NPATR? 

XW(J,2)= XW(J,?> / 0 
YW(J,2)= YW(J,2) / D 
CURV(J»2)=CURV(J.2)*0 
CONTINUE 

Rhooi=ai*poi/rg/toi 

UREFP*SQRT( GC*RG*T01) 

K.jft is NUHRFR of POINTS CLOSE TO .lET NOZZl E WAl I S 
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$NA1300 

SNA1310 

5NA1320 

5NA1330 

5NA1340 

5NA1350 

SNA1360 

SNA1370 

5NA1380 

SNA 1390 

SNA1400 

SNA1410 

INA1420 

SNA1430 

SNA1440 

JNA1450 

SNA1460 

5NA1470 

$NA1480 

SNA1490 

SNA1500 

SNA1510 

SNA1520 

$NA1530 

SNA1540 

INA1550 

SNA1560 

$NA1570 

SNA1S80 

SNA159C 

SNA1600 

SNA 16.10 

SNAlftao 

SNA1630 

SNA1640 

fNA16S0 

SNAl 660 

SNA1670 

5.NAI680 

SNA1690 

SNA2700 

SNAi?lO 

SNA1720 

5NA1730 

SNA1740 

$NA1750 

SNA1760 

SNA1770 

SNA1780 

SNA1790 

SNAieOO 

SNA1810 

SNA1820 

SNA1830 

$NA1840 

SNA1850 

$NA1860 

SNA1870 

$NA1880 

SNA1890 

SNA1900 

SNA1910 

5NA1920 

SNA1930 

SNA] 940 
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L 


« a 


kwall is number of points close to 

KJET= 6 
KWALL= 20 

N= k 3 ♦ Kl + k2 ♦ 2 * kwall ♦ A * 
NN=N-1 
NX1 = 1 

NX2=K3*KWALL*KJET+1 

NX3=NX2 

NX4=N-K2-KWaLL-KJeT 
NX5=NXA 
NX6=N 
NX (1)=NX1 
NX<2)=NX2 
NX(3)=NX3 
NX(A)=NX4 
NX(5)=NX5 
NX (6)sNX6 
Ll= NX(3) 

L2= NX *4) 

NMAX= (LI ♦ L?) / 2 
NF= 1 
P03=P02 
T03=T02 

massi=ai*massi/rhooi/upefp/d 

MASS?=A1*MASS2/RhO01/UREFP/O 


DUCT WALLS 

SNA1950 


SNA1960 


SNA1970 


$NA1980 

KJET ♦ 1 

5NA1990 


SNA2000 


$NA2010 

SNA2020 

5NA2030 

5NA2040 

SNA2050 

SNA2060 

$NA2070 

$NA2080 

SNA2090 

5NA2100 

SNA2110 

SNA2120 

SNA2130 

SNA2140 

$NA2150 

SNA2160 

$NA?170 

$NA2180 

SNA2190 

$NA2200 
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e 

e 

e 


28 


29 

921 

922 

923 

924 

925 

926 

927 


SNA2210 

SUBROUTINE INCOMD IS USED TO CALCULATE THE INITIAL CONDITIONS 5NA2220 

SNA2230 

CALL lNcONn<N*NN*Kl »K2fK3»KWALL*KJET*Ll*L2t $Na22^0 

1 MaSS1*MaSS2*AMASS1*AMaSS2*AMASS3» SNa2250 

2 RG*GC*G*XN0Z*YN0Z*TN0Z»XW*YW*NPA1R1,NPAIR2) $NA2260 

AMA1=AMA5S1*UREFp*O*RHO01/AI SNA2270 

AMA2=AMASS2*UREFP*D*RH001/A1 5NA2280 

AMA3=AMASS3*UREFp*D*RH001/A1 SNA2290 

WRITE(6»921) SNA2300 

62=0/ (G-1.) SNA2301 

VNOZI=UftEFP*SQRT(2.*62*(l.-<P02/P01>**a*/G2) > ) SNA2302 

TN0ZI=AMA1*VN0zI SNA2303 

WRITE (6*927) TNOZI SNA2304 

WRITE (6*922) AMAl * AMA2*AMA3 INA2310 

WRITF(6,923) SNA2320 

Do 28 J=1*N SNA2330 


XS=X ( J)*0 
DNS=ON( J)*D 
YS=Y( J)*D 

TX=THETA(J,1 ) *57,2958 
XR=R( Jf n/D 

write ( 6*924) J,XS,YS*DNS*TX,XRtPSl (J> , 

WRITE(6*925) 

DO 29 J=1 *N ' 

POD=PO ( J)*P01 

PSD=PS( J, l)*POl 1 

T0=TS(J*1)*T01 

KHOD=RHO ( J* 1 ) *PHOOI 

UV=U( J*I)*UREFP 

WRITE<6, 26) J,PnD*P50*TO*RriOD*UV*M( J) 

FORMAT (IHl •///*2^X,18HINITIAL CONDITIONS//) 

format (7X *BHAMaS51 =F I 0 *5 *5X , 8HAMA$S2 =F 1 0 .5* 5X , RHAMASS3 =F10,5/) 
FORMAT (7X ,1hJ*9X* lHX*9X*lHY*ex*2HDN*6X*SHTHETA*8X*lHK*7X*3H0MG/> 
Format (5x, 13 *2x*4F 10.3 *F1U.6*F10 *4) 

format (///7X*1HJ,8X*2HP0*8X*2HPS*8X*2HTS*7X*3HRH0*9X*IHU*9X*1HM/) 

FORMAT(5x*I3*2x,?£l0.3,F10.1*Fl0.5,F10.1*Fl0.3) 

format (7X*41HISEvTP0PIC nozzle thrust per UNIT WIDTH =*ElO*3) 

DO 30 J=1*N 
theta ( J*?')=THETA( J*l) 

PS( J,?)=PS( J*1 > 

TS( J*2)=TS(J*I ) 

RH0(J,2) =RH0(J*1> 

U( J*2)=U( J*I) 

R( J,?)=R(J*1) 


5NA2340 

$NA2350 

$NA2360 

1.NA2370 

5NA2380 

SNA2390 

INA2400 

iNA2410 

$NA2420 

SNA2430 

SNA2440 

$NA2450 

4NA2460 

$NA2470 

$NA2480 

5NA2490 

SNA2500 

SNA25I0 

5NA2520 

5NA2530 

SNA2531 

$NA2540 

$NA2550 

SNA2560 

$NA2570 

SNA2580 

$NA2590 

SNA2600 
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30 C(J) s 0. 

DEL1=0.0 

DEL2=0*0 

DEL3=0. 

DEL4=0. 

COrE3sDN (NX2) -on (NXI ) 

C0rE1=DN(NX4)-DN(NX3) 

C0rE2=0N (NX6) -ON (NX5) 

DPDSAs -0.00 1 
USTARI = 0.0 
USTArO s 0.0 
UGOsO. 

OGisO.O 
NGRIDsO 
NPRINT= 0 
DSAs .02 * ON(N) 

CC* .07 
NC0RE*1 
NSEPsO 
R0ECAY=4. 

ICOUNT^O 

c 

c 

c 

C BEGINNING OF MARCHING CALCULATION 

C 

c 

c 

00 4500 NSTEP = ItlOO 
ICHS = 0 

c 

c find middle of Duct 

c 

DO 527 J“l*N 

IF (DNIJ) .LT. DN(N)*0.5) N0UCT= J 
527 CONTINUE 

IF INSTEP .EQ. 1> GO TO 77 
C 

C SUBROUTINE EDDY CALCULATES THE EDDY VISCOSITY FOR SUB. SOLV 

C 

Call EODY(N,NN»Nx*U.PS.CC,DS»E*EI,RHO,viS*R*DNt 

I S4, S5» 0EL1*DEL2»0EL3,DEL4,C0RE1,C0RE2.I:0H£3) 

77 CONTINUE 
El(l» = 0. 

00 40 J=2*NN 
El(J) » RhOCJ»1)*U<J»1) 

AO Continue 
E l <l^)sO. 

c 

C move to next point on WALL 

c 

OS(NMAX):s OSA * (1.04) *» (NSTEP-D 

IF (OS(NMAX) .lT, ,02 ♦ Dn(N) ) DSiNMAX)* ,02 * DN(N) 

c 

c Calculate os(I) values 

C (MIDDLE OS value CALCULATED FIRST 

C then os values calculated out TO BOTH WALLS 

c 

NPRRsNMAX^I 
DO 50 J=NPRR*N 
JMs J-1 

C1*RH0(J*2) ♦RHO( JM*2) 

C2=U(J*2)^U(JM,2) 

C3=R(J*2)+R(JMf2) 

C4= PSI (J) - PSI (J-l) 
fs * rA*c:i/(Ci*r2) 


SNA2610 

SNA2620 

SNA2630 

SNA2640 

SNA2650 

$NA2660 

SNA2670 

SNA2680 

SNA2690 

SNA2700 

SNA2710 

SNA2720 

SNA2730 

5NA2740 

$NA2750 

5NA2760 

SNA2770 

SNA2780 

5NA2790 

$NA2800 

SNA2810 

$NA2820 

SNA2830 

5NA2B40 

$NA2850 

SNA2860 

SNA2870 

SNA2880 

SNA2890 

$NA2900 

SNA2910 

$NA2920 

SNA2930 

$NA2940 

SNA2950 

$NA2960 

SNA2970 

$NA2980 

5NA2990 

SNA3000 

SNA3010 

$NA3020 

$NA3030 

SNA3040 

INA3050 

SNA3060 

SNA3070 

$NA3080 

SNA3090 

$NA3100 

$NA3110 

SNA3120 

5NA3130 

5NA3140 

SNA31S0 

$NA3160 

5NA3170 

$NA3180 

$NA3190 

SNA3200 

SNA3210 

SNA3220 

SNA3230 

SNA3240 

SNA3250 

SNA3260 
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50 


52 

C 

C 

C 

C 


C 

C 

C 


C 

c 

c 


400 

C 

C 

C 


51 


C 

C 

C 

C 

C 

C 

c 

c 

c 


c 

c 

c 


DS(J) s (U+C5)/(l.-C5>*DS(JM> 

continue 


nprr=nmax-i 
DO 52 MM=lfNPRR 
J= nMAX - HM 
JM= J*1 

Cl=RHO(Jt2) ♦RHO(JM#2) 

C2*U( J,2) ♦U(JH*2) 

C3» R(J*2) ♦ R<JM*2> 

C4* PSI(J) - PSI(JM> 

CS s C4*C3/(CI*C2> 

DS(J> <1.*C5)/(1.-C5>*DS(JM) 

CONTINUE 


$NA3270 

$NA3280 

JNA3290 

INA3300 

5NA3310 

SNA3320 

5NA3330 

SNA3340 

SNA3350 

$NA3360 

SNA3370 

SNA3380 

$NA3390 


- SNA3400 


HOVE TO NEW POINT X(l) AND X<N) ALONG WALL SURFACES SNA3410 

AT A DISTANCE SPECIFIED BY DS VALUES JNA3420 

$NA3430 

0X1= X(l> INA3440 

CALL OSMOVE(XW«YW»xn)tY(l)«OS(l)*THETA(l»l)» R(l*l)t NPAIRlfl) $NA3450 
CALL OSMOVE(XW«YW«X(N) tY(N) « DS(N)« THETA(N*l>t R(N*I)t NPAIR2*2> SNA3460 
1F(X(1) .GT. XWiNPAlRl •!)) STOP SNA3470 

IF (X(N> .GT, XW(NPAIR2 ,2) ) STOP SNA3460 

0X1= X(l) - DXl SNA3490 

SNA3500 

COMPUTE Y* theta* AND R VALUES CORRESPONDING TO nEW X VALUES 5NA3510 


Call walls(xw*yw* xdi* Y(D* thetai* Rit npairi* d 
Call walls(xw*yw* x(N)» y<n)* thetaz? r2» npaIr2« 21 

Y2= Y(N) 

Compute curvature R(j*2> 

DO 400 J=1«N 
Owl= DN(J) 

DW2= DN(N) - DN(J) 

K( J*2)=R1*U.-DW1/0N(N) )*EXP(-ROECAY*DWl*ABS(Rl) >*R2*0W1/0N(N>* 
1 EXP (-RDECAY*DW2*ABS(R2) > 

CONTINUE 

COMPUTE THETA (J*2) 

Oo 51 J=2*NN 

THETA(J»2)= THETAI ♦ (THETA2-THETA1 ) * ON(J)/On(N) 

Continue 

ThETA{1»2)= THETAI 
ThETA(N*2)= THETA2 


THIS SECTION attempts TO SATISFY CONTINUITY 
LOOKS FOR A PS<1*2) SUCH THAT Y2 - Y(N) * 0«0 
DPDSA* DPDSB are PRESSURE GRADIENTS 


IF(A8S(0PDSA) .LT. l.E-08) OPDSAs>*001 
DPDSB= DPDSA * 0.9 
DDISTB= 1.0 
DO 650 ITER= 1*40 

IF (ITER .EQ. 1) PS(1*2>= PS(1»1) DPDSA 

AT THIRD ITERATION GUESS NEW DPDSB 

USING EXTRAPOLATION THROUGH OPDSA AND PREVIOUS DPDSB 
IF (ITER *GE. 3) 

1 DPDSBa (DDISTB*OPDSA - 0DISTA*0PDS8) / (DDISTB-DDISTA) 
IF (ITER .GF. ps(i.p)= psn*n ♦ nPUSB 


SNA3520 
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SNA3540 

INA3550 
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SNA3640 
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$NA3660 

$NA3670 

SNA3680 

$NA3690 
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5NA3720 

5NA3730 

$NA3740 

SNA3750 

$NA3760 

SNA3770 

$NA3780 

INA3790 
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SNA3810 
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INA3830 
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SNA3850 
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$NA3870 

SNA3880 

SNA3890 

SNA3900 
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DO 60 J*?*N 

SNA3930 


JM= J-1 

SNA39A0 


C?=U(J*2)->U(JMf2) 

SNA3950 


C3=R<J»2)+R(JV*21 

SNA3960 


C4s PSKJ) - PSI(J-n 

SNA3970 


PS(Jf?)=PS(JM,2)^ C2«C3*C4/4.0 

$NA3980 

60 

CONTINUE 

SNA3990 

C 


SNAAOOO 

C 

SUBROUTINE SOLV IS USED TO SOLVE UtT ON M«2 LINE 

SNAAOlO 

C 


SNAA020 


call SOLV<DS*N*NNfPRltPRTlfG) 

SNAA030 


U(l«2) > 0* 

SNAAOAO 


DO 70 J * 2*NN 

SNAAOSO 


U(J»2) = H(J-l) 

SNAA060 

70 

CONTINUE 

SNA4070 


U(N«2) = 0. 

SNAAOeO 


call S0LV(0S*N*NN#PR*PRT*G) 

$NAA090 


00 80 J = 2*NN 

SNAAIOO 


TS(J*2) = H(J-l) 

SNAAIIO 

BO 

CONTINUE 

SNAA120 


TSdf?) = TS(?*2) 

SNAA130 


TS(N*2) = TS(N-1*2) 

SNAAIAO 


DO 90 I = 1*N 

SNAA150 


RH0<I*2) s PS(I*2)/TS(I*2) 

SNAA160 


. RHO(Itl) sRH0(I>2) 

SNA4170 


VIS<I)*<TREF^A2)/<T01*T5(I*2)*A2) 

SNAA180 


VIS(I)=VISm*(T0l*TS(If2)/TREF)**0*5 

SNAA190 


V I S ( 1 ) sV I S < l i *RH0REF*NUREF/ <P0 1 *S0RT < GC/RG/TOl > ) 

SNAA200 

90 

CONTINUE 

SNAA210 

C 


SNAA220 

C 

GIVEN PS (1*2)* CALCULATE UPPER WALL COORDINATE Y(N) 

SNAA230 

C 


SNAA2A0 


Call cords (N*XfY*THETA»U*RHO*P.SI»ONtNDUCTfDOlST) 

SNAA250 

c 


$NA4260 

C 

IF Y2 .EQ. VCNl THEN THIS PS<1*2> SATISFIES CONTINUITY 

SNA4270 

C 


SNAA280 


IF (ITER .EO, X) ODISTA* DOIST 

5NAA290 


IF (ITER ,GE. 2) 00IST8= OOIST 

SNAA300 

C 


SNAA310 

C 

IF SUFFICIENTLY SMALL INTERVAL? THEN EXIT 

$NA4320 

C 


SNAA330 


IF (ABS(DOISTB) .LE. O.l ** A) GO TO 660 

SNAA3A0 

650 

continue 

SNAA350 


WRITE (6t916) OPOSA*OPDSBf ODISTBfNSTEP 

SNAA360 

916 

Format (/7X»23HN0 convergence? DPDSA *?F10.5?BH DPDSB b,F10«5? 

SNAA370 


1 9h DDISTR »?F10*5?6H ST£P=I4) 

SNAA380 

660 

continue 

SNAA390 

C 


SNAAAOO 

C 

SAVE PRESSURE GRADIENT THIS STEP 

SNAAAIO 

C 


SNAAA20 


OPDSA= DPDSB 

5NAAA30 


GisG-J. 

SNAAAAO 


Gp=G/Gl 

SNAAA50 

C 


SNAAA60 

C 

COMPUTE MACH NUMBERS M(I) 

SNA4A70 

c 


SNA4A80 


DO 220 I=l?N 

SNAAA90 


M(I)=U(I,2)/S0RT(G*TS(I?2)) 

SNAA500 


GA = 1, ♦ Gl/2.*(Mn)*M(IM 

SNAA510 


P0(I) - PS(I?2)*GA**G2 

SNAA520 

220 

Continue 

SNAA530 

C 


SNAA5A0 

SNAA550 

C 

c 

c 

REMOVE EXTRA GRID POINTS? FIND EDGES OF FLOWS 

SNAA560 

SNAA570 
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Call REM0VEINtNN*NGRIU«NX«NMAX«PSitSl»S2»S3*S4«S5*D!>*VIS*DN* 

1 X*Y*MfPO*H*7HETA*K»PS«U*TS*RHO) 

CALL LOOK (N*NN,K3tNC0RE*DELl *0EL2,DEL3tDEL4*C0REl *C0RE2*C0RE3» 

1 NX*PO«DNfU«D1P«NMAX> 

ShRIN=U(?,2)*VIS<2)/<DN(2)-DN(1) ) - ,5* (PS ( 1 *2) -PS { 1 * 1 ) ) /OS ( 1 > * 

1 (DN(2>-0N( 1) ) 

SHR0R=U(NNt2>*VfS(NN)/(DN(N)-DN(NN) ) -.5*(PS(N»2)-PS(Ntl ) »/ 

1 0S(N)*(DM(N)>DN(NN) > 

IF( (SHRIN.LE.O.) .OR. (SHROR.LE.O.) ) NSEP=1 
1F(NSEP.E(3.1> NPRlNTsNPCYCL-1 
IF(NSEP.EQ.I) ICOUNTsl 
C 

C — 

C SAVE VALUES COMPUTED THIS STEP 

C — 

C 

Do 222 1=1 fN 

PS(I*1) « PS(1«2> 

U(I.1> =U(I«2) 

TSd.l) * TS(I,2» 

ThETA(I*1)=THETA(I»2) 

222 R(Ifl) = RU«2) 

C 

C 

C OUTPUT section 

C 

C 

DSlOs DS(1) * 0 
DS20= DS(N) * 0 
DSUM1= DSUMl ♦Oslo 
OSUM2= OS2Q ♦ DSUM2 
NPRINTs NPRINT*! 

IF ( (X(1)*D) ,GE# XPROF(nF) ) NPRlNT= NPCYCL 

lF(X(l)*D.GE.XPROF(NF) ) ICOUNT=l 

IF { <X(1)*D) .GE. XPROF(nF) ) NF* NF ♦ I 

IF (NSTEP .£0. 1) NPRINTs NPCYCL 

IF (NPRINT .LT. NPCYCL) 60 TO 45 

NPRINT= 0 

0X0= 0X1 * D 

X 1= X(l) *0 

X0= X(N) * 0 

PH20I=(PSn.2)*POl-P02)/A3 
PH200= (PS(N*2)*P0l-P02)/A3 
IF (SHRIN .GT. 0.0) 

1 USTARI=SQRT(SHRIN/RH0(2*2) )*UReFP 

If (SHROR .GT. 0.0) 

1 USTAR0=SQRT(SHR0R/RH0(NN.2) )*UREFP 

Hl= H(1.2) / D 
R0= R(N.2) / 0 

RnI= 2.0 * (R(2.1) / RH0(2.D) 

RNI=RNI/(S5(2)*(U( 3*1)-U(2. 1) )^S4(2)*(U(2*1)-U("T*1) ) ) 

RnO= 2.0 * (R(NN*1) / KHO(NN.l)) 

RNO=RNO/ < S5 ( NN ) * ( 1) ( N * 1 ) ->U ( NN « 1 ) ) +S4 ( NN ) * ( U ( NN f 1 ) -U ( NN- 1 • 1 ) ) ) 

UMAX= U(NMAX*2) • UREFP 

XTSI= TS(1.2) * TOl ♦ 0.S/G2 * TOl * (U(1.2>)*<*2 
XTSO= TS(N.2) * TOl ♦ 0.S/G2 * TOl * <U(N*2))**2 

INRITE (6*908) NSTFP*XI*PH2oI*USTARI.RI*RNI.XO*PH200.USTARO*RO*RNO* 
1 UMAX. ITER 

906 F0RMAT(lMl«5X.5HNSTEP.6X«2hXI«8X*2HPI.4X*6HUSTARr*8X*2HKI*7X* 

1 3MRNI t8X*2HX0.8X*2HP0*4X,6HUSTAR0.8X*2HK0*7X*3HRN0*6X*4HUMAX/ 

2 5X. 13*2 (F10.4.EI1.3.2F1 0.4. FIO.S) *F10.1//7X« 

3 22HNUM8ER OF ITERATIONS =*l3/> 

WRITE (6*907) 0X0. DS1Q.DS2Q.0SUMI,DSUM2 

907 F0RMAT( 7X* 9HD£LTA X =*F10.S*12H DELTA Sl =*Fl0.5. 

L 12H delta SO =*F10.5*6H SI =F10.5*6H SO =*F10.5/) 

Sho TMbSHR TN»PO 1 


INA4590 
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5NA4620 

SNA4630 

SNA4640 

$NA4650 

SNA4660 

SNAA670 

SNA4680 

SNA4690 

SNAA700 

SNA4710 

5NA4720 

SNA4730 

$NA4740 

SNA4750 

1NA4760 

SNA4770 

$NA4780 

5NA4790 

$NA4800 

$NA4810 

SNA4820 

5NAA830 

SNA4840 

SNA4850 

$NA4860 

SNAA870 

$NA4880 

fNAA890 

$NA4900 

5NA4910 

$NA4920 

SNA4930 

£NA4940 

$NA4950 

$NAA960 

$NA4970 

SNA4980 

JNA4990 

SNASOOO 

SNA5010 

$NA5020 

INAS030 

iNA5040 

SNA5050 

SNA5060 

$NA5070 

SNA5080 

SNA5090 

SNASlOO 

SNA5110 

$NA5120 

$NA5130 

$NA5140 

SNA51S0 

SNA5160 

5NA5170 

JNA5180 

SNA5190 

$NA5200 

$NA5210 

SNA5220 

$NA5230 

SNA5240 
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SHROR=SH«OR*POl 

WRITE (6tQ12) SHRiNfSHROR 

TMOMXsO. 

NNY=N-1 
00 95 J=1*NNY 

95 Tm0MX=TM0MX+(PH0( J»1)*U(J*1)**2*C0S(THETA(J,1) ) ♦rHo(J*1 *D* 

I U( J*lf 1)**?*C05(THETA(J+1*1) ) )/2.*(Y(J*U»Y(J) ) 

tmomx=tmomx*poi*d 
CT=TM0MX/TN07I 
WRITE<6»931) TMOMX*CT 

912 FORMAT! 7X*7HSHRlN =E1 2. 5 * 5X , 7HSHR0R =El2.5/> 

931 FORMAT (/7X»37hTOTA|. AXIAL MOMENTUM PER UNIT WIDTH =fEl0.3//7X* 

1 iahthrust ratio =*E10.3) 

IF (NQUICK .GT, 0> GO TO ^^5 
IF (ICOUMT.EQ.O) GO TO 
ICOUNT=0 

C0RE1Q= COREl * D 
CORE?Q= C0RE2 * 0 
C0RE3Q= C0RE3 * 0 
DEl 10= OELl * D 
0 El2Q= DEL2 * D 
0EL3Q= 0EL3 * 0 
0EL40= DEL4 * 0 

WRlt£(6»909) (NjX(I) tl=l*6) »C0RE:Q»C0RE2Q*C0RE3 Q,DEl1Q»DEl2Q* 

I DEL3QfDELAQ 

9u9 format ( 7x*2HNl,5x*2HN2f5x»2HN3,5X*2HN4,5Xf2HN5,5X»2HN6/5x* 14, 

1 5(3X,14)// 7X,5HC0RE1 ,5X,5HCORE2*5X,5HCORE3*/4X,3F10,5// 

2 7X,4HDEL1 *6X,4HDEL2,6X,4HOEL3,6X*4HDEL4/4X,4F10.5//) 
WRITE <6,910) 

910 Format ( 7x,l Hj,7X,3hl/R, 5X»5HTHET A, 9X,1HX,9X,1HY,6X,4HYREL>6X, 

1 4HURFL,5X,5HdOREL,8X,2HPS,9X*1hE,5X,6HTOTEMP/) 

00 100 J=1,N 

THEI0=THETA( J»2)*1B0./3.1^16 

XR=R(J,2)/D 

YS=Y<J)*D 

YREL= Y(J> / Y(N) 

XS= X(J)*D 

UV=U(J,2) / U(NMaX,2) 

POO=PO ( J) 

PSD=i^S ( J,2) *P01 

To= TS(J,2) » Tol ♦ 0,5/62 * TOI * (U(Jt2))**2 

TOMAX= TS(NMAX»?) * TOl ♦ 0.5/G2 * TOI * (U (NMAX, 2) ) **2 

XTS= (To - XTSO) / (TOMAX - XTSO) 

XT=E < J) *UREFP*D/S0RT (Al ) 

OND=DN ( J) *D 

write ( 6,911) J,xR,THEIO,XS*YS,YREL,ijV,POO,PsD,XTf TO 

911 format ( 5X, 13, 2X, 1^1 0.5, 3F1U. 4, 3F 10.5, El 1,3, FI 0*5, FI 0,2) 

loo CONTINUE 

45 CONTINUE 

IF (NSEP.EO.l) GO TO 300 

4499 Continue 

4500 Continue 

STOP 

300 WRITF(6,930) 

930 F0RMAT(7X,24(1h*)/7X,38HNEGATIVE SHEAR STRESS AT THIS STATION./ 

1 7X,36H THIS indicates POSSIrlE SEPARATION./ 

2 7X,48H subsequent RESULTS SHOULD BE USED WITH CAUTION,// 

3 7X,24(1H*)> 

NSEP«0 

GO TO 4499 
END 


$NA5250 

$NA5260 

$NA5261 

SNA5262 

SNA5263 

SNA5264 

$NA5265 

5NA5266 

$NA5267 

SNA5268 

$NA5270 

SNA5271 

SNA5272 

SNA5280 

SNA5290 

5NA5300 

$NA5310 

SNA5320 

SNA5330 

5NA5340 

$NA5350 

$NA5360 

SNA5370 

SNA5380 

$NA5390 

SNA5400 

$NJA5410 

SNA5420 

SNA5430 

$NIA5440 

SNA5450 

SNA5460 

SNA5470 

$NA5480 

SNA5490 

SNA5500 

$NA5510 

$NA5520 

$NA5530 

$NA5540 

$NA5550 

SNA5560 

$NA5570 

$NA5580 

$NA5590 

5NA5600 

SNA5610 

$NA5620 

SNA5630 

SNA5640 

SMA5650 

$NA5660 

$NA5670 

SNA5680 

SNA5690 

INA5700 

SNA5710 

$NA5720 

$NA5730 

SNA5740 

SNA5750 
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ooon ooonoono oo o or>oor>no r>..o o o o o 


c;iJBPOUTlME INC0N;')(M»NtvJ»Kl,i'v£;,K3»KW/\LL*KJL» »LJ »Lif, 

1 MASSl ♦MASSa * amass 1 tAMASS2»AMASS3* 

2 PG*r,C»GtXNOZ,YN0Z»TNO2»XW*YW*NPAlRl,NPAlR2) 

SlIRROUTINE INCONO 

Calculates initial conditions 
Input variables 

NPAIRi: NUMBER OF POINTS ON LOWER WALL 

NPAIR?: number of points ON UPPER WALL 

XW(I*D* Ywn*l>: data points specifying LOWER WALL 

xw(i»?>* YW(i«2): data points specifying upper wall 

XI : starting point on lower wall 

XC»YC: centre of initial RADIUS OF CURVATURE 


REAL MASSl f MASS?*NU0*M (1901 

dimension XW(99. 2) *YW(99»2) *SK190) ,52(190) ♦S3U90) * S4(190)» 

1 S5(190) .E(190) ,E1 (190) ,DS<190) ,X(190> ,Y(190) ♦ 

2 P0(190) ,H(190) , THETA (190,2) ,R ( 190,2) ,U ( 190 ,2) , 

3 TS(190,?) ,RHO(190,2) ,VIS(190) 
dimension OMG(I90) »ON(190) ,PS(190,2) 

OniiBLF OMG(190) ,0N ( 190 ) ,PS ( 190 ,2) 

DIMENSION PSI ( laO) 

DOUBLE PSI (IPO) 

Double xt,p3,py,pp,dy,psi , nPsi,ZY,DsQRT,0Ps,z 
rOMMON/SOLV/ SI ,S2,S3,E,E1 ,H,SA,S5 
COMMON/INEDY/PSI »RHO,VIS 

COMMON/INCO/ U»PB,TS,M,PO,THETA,R,x,Y,DN 

COMMON /CONST/ P0l,P02,PC3, T01,T02,T03* XI, XC, YC 

DIST (XX,YY,XXC,YYC) =SORT ( ( ABS ( XX-XXC ) *<»2) ♦ OABS ( YY-YYC ) **2) ) 

LOOK FOR ARC TANGENT TO BOTH WALLS WHICH PASSES THROUGH NOZZLE 
XI AND X? ARE STARTING POINTS ON LOWER AND UPPER DUCT WALL 

XNOZ. YNOZ IS middle OF NOZZLE 
FIND COORDS OF EACH EDGE 

nozzle width is 1.0 

XN0Z1= XNOZ ♦ { (SIN(TNOZ)) * 0-5) 

XN072* XNOZ - ( (SIN(TNOZ)) * 0.5) 

YN0Z1= YNOZ - (flBS (COS (T mOZ) ) * 0.5) 

YMOZ?= YNOZ ♦ ( ABS (COS (TmOZ) ) * 0.5) 

COMPUTE ARC LENGTH ACROSS DUCT (ON?) 

ARC lengths to nozzle (DNI, DNO) 

CALL ARCDIS (XN021 ♦ YNOZ 1 , TnOZ , X 1 , RC 1 , X W , YW ,MPA I R 1 , 1 ) 

CALL WALLS (XW,YW,X1 ,Y1 ,T1 ,R1 ,NPAIR1,1) 

IF (PCI .LT. 100,0) DNI= RCl * ARS(TNOZ-Tl) 

IF (RCl .GE. 100,0) DNI= DIST (XNOZl , YNOZl ,Xl , Y1 ) 

CALL ARCDIS(XMD22,YN0Z2,Tn0Z,X2,RC2,XW,YW,NPAIR2,2) 

CALL WALLS(XW,YW,X2,Y2,T2,R2,NPAIR?,2) 

IF (RC2 .LT. 130.0) DNO= RC2 * AbS(TN0Z-T2> 

IF (PC2 .GE. 100. 0) ONO* DIST (XN0Z2, YN0Z2,X2, Y2) 

DN2= DNI ♦ DNO ♦l.O 

X(1)=X1 

Y(l)=Yl 

X(N)= X2 

Y(N)=: Y? 

THETA(1,1) s T1 
THETAIN*!) s T? 


SIN 

0 

SIN 

10 

SIN 

20 

SIN 

30 

SIN 

40 

SIN 

50 

SIN 

60 

SIN 

70 

SIN 

80 

SIN 

90 

SIN 

100 

SIN 

no 

SIN 

120 

SIN 

130 

SIN 

140 

SIN 

150 

SIN 

160 

SIN 

170 

SIN 

180 

SIN 

190 

SIN 

200 

SIN 

210 

SIN 

220 

SIN 

230 

SIN 

240 

SIN 

250 

SIN 

260 

SIN 

270 

SIN 

280 

SIN 

290 

SIN 

300 

SIN 

310 

SIN 

320 

SIN 

330 

SIN 

340 

SIN 

350 

SIN 

360 

SIN 

370 

SIN 

380 

SIN 

390 

SIN 

400 

SIN 

410 

SIN 

420 

SIN 

430 

SIN 

440 

SIN 

450 

SIN 

460 

SIN 

470 

SIN 

480 

SIN 

490 

SIN 

500 

SIN 

510 

SIN 

520 

SIN 

530 

SIN 

540 

SIN 

550 

SIN 

560 

SIN 

570 

SIN 

580 

SIN 

590 

SIN 

600 

SIN 

610 

SIN 

620 
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c 

c 

c 

c 

c 

c 

c 

c 


c 

c 

c 

555 


C 

C 

C 


c 

c 

c 

c 


SOLVe FOP MASS SPLIT {N027LE EXIT CORRESPONDS TO ON(Ll) ) 

DUCT COMPUTFO MUST CORRESPOND TO ACTUAL DUCT WIDTH 
THUS DNU) SHOULD EOilAL 0*0 

FIND LOCATTO^f OF NOZZLE 

(DIFFERENT CALCULATIONS FOR SPlMAX DEPENDING ON 
NOZZLE LOCATION) 

TF (DNI ,LT. DNO) NZONE= 0 
IF (DNI ,GE. ONO> NZONEs 1 
IF (NZONF .EO, 0) NMIO= L? 

TF (NZONF *EO. 1) NMID= LI 
Gls 1.0 / (G-1.0) 

G?= G/( G-1.0) 

PDECAY=A. 

A=DNI* I .♦DNO 

IF(NZONE,EO.O> RTEST=R1*( J .- (DNI ♦ I . ) /A ) «EXP (-RDECAY* (DNI ♦! • ) * 

1 APS(R1 ) )^R?*(nNI+I,)/A*EXP(-RDECAY*DNO*ABS(DNO)) 

IF(NZONE.EO.l) RTEST=Rl*n .-DNI/A)*EXP(-RDECAY*D nI*ABS(R1) )♦ 

1 R2*DNI/A*EXP(-RnECAY*(nN0^1.)*ARS(R2) ) 

TMAX= (T02/T01 - .0001*( (G-1.0)/(G*2.0) ) ) 

1 * (P03/P01) ** ((G-l.f!)/G> * (T01/T02) 

»MNMID>= SQRT(2.0 * G1 * (1.0/ TmaX - l.O) ) 

oSNMIO= TMAX ** G2 

PH0(MMID,1)= POP / POl * TOl / T02 

H2MIN=ABS(.0001-?.*G*PSNMID*(M(NM1D) )**2*RTEST/RH0(NMID*1) ) 
M2MIM= SQRT(U?MTN) 

SPLMAX= l.O - ( P02 / POl * TOl / T02 * U2MIN * DNO / MASS2) 
SPLMAX= ABS(SPLMAX) 

IF (SPLMAX .CT. 0.99) SPlMAX* 0,99 

THIS IS ENTRY POINT FOR LOWERING SPLMAX IF TMAX iS EXCEEDED 
CONTINUE 

SPLMIN= 1.0 - SPLMAX 
SPLTST* DNI / (ONI ♦ DNO) 

XA= SPLTST 

IF ((NZONE .FO. 0).AND.(XA .LT. SPlMIN*! • 1 1 ) ) XAs SPLMIN*1,11 
IF ((NZONE .EO. D,AN0.(XA .GT, SPlMAX* 1,0)) XA* SPLMAX* 1.0 
XB= XA * 0.9 
niFFB= 1.0 

NMID= (LI ♦ L2) / 2 

DO AO ITFR= 1.50 

IF (ITER .EO. 1) SPLIT=XA 

IF (ITFR .GE. 3) XB* (DlFFB*XA - DIFFA*XB) / (DIFFB-DIFFA) 

TF (ITER .GE, ?) SPLIT* XB 
AMASSlsMASSl 
AMASS3* SPLIT * MASS2 
AMASS?* (I.O-SPLIT) * MA$S2 

stream FUNCTIONS 

OMGN= AMASSl ♦ AMASS2 ♦ AMASS3 

0M6J3S AMA5S3 

OMGJ2* AMASS3 ♦ AMASSl 

OMJET IS THE SMALL DELTA oMG AT JET WALLS 

DMWALL is THE SMALL DELTA OMG STARTING IN FROM DUCT WALLS 

DMWALL* .00001 * OMGN 
DMJET* DMWALL * 100.0 
DOMG3* 0MGJ3 / FLOAT (K3) 

O0MG2* (OMGN - OMGJ?) / FLOAT (K2) 


SIN 630 
SIN 640 
SIN 6S0 
SlN 660 
SIN 670 
SIN 680 
SIN 690 
SIN 700 
SIN 710 
SIN 720 
SIN 730 
SIN 740 
SIN 750 
SIN 760 
SIN 770 
SIN 780 
SIN 790 
SIN 800 
SIN 810 
SIN 820 
SIN 830 
SIN 840 
SIN 850 
SIN 860 
SIN 870 
SIN 880 
SIN 890 
SIN 900 
SIN 910 
SIN 920 
SIN 930 
SIN 940 
SIN 950 
SIN 960 
SIN 970 
SIN 980 
SIN 990 
SINIOOO 
SINlOlO 
SIN1020 
SIN1030 
SIN1040 
SIN1050 
SIN1060 
SIN1070 
S1N1080 
S1N1090 
SINl 100 
SINIIIO 
SIN1120 
SIN1130 
SIN1140 
SIN1150 
S1N1160 
SIN1170 
S1N1180 
SIN1190 
SIN1200 
SIN1210 
SIN1220 
SIN1230 
S1N1240 
SIN1250 
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c 

c 

c 

<►0 

42 

902 


OnMGl X 
psi n ) = 
PSI (Ll 
PSI (L2>= 
PSI <N)= 
CALL 
CALL 
CALL 
CALL 


(OMGJ? - 
0.0 
0MGJ3 
OMGJ? 
OMGN 


OMGJ3) / FLOAT (Kl> 


IF 

IF 

IF 

IF 


OMGSET (PSI, 1 .LI »KWALL.KJET,0MWALL*DHJFT*D0MG3.RST3*RFTN3) 
0MGSET(PSI.L1*L2»KJeT.KJET»DMJET.DMJET»D 0MG1*RST1.RFIN1> 
OMGSET ( PSI »L?»N.K JET »KWALL. DM JET ♦0MWALL*D0MG2.RST2.RFIN2) 
TMPSFT (N.NN»RG»GC.G»0N2*R1*R2.L1 tL2*0NI*DN0* 

PSIt R. M« TS» PS« POt RHOt U« ONvMAXERR) 
((MAXERR .EO. 1 ). and. (ITER .EQ. 1)) GO TO 589 
((MAXERR .EO. 1). AND. (ITER .NE. 1)) GO TO 599 
(ITER .EQ. 1) DIFFAs DN(1) 

(ITER .GE. ?) DIFFRs DN(U 


CHECK FOR termination 

IF (ABS(DIFFB) .L£. O.OOl) GO TO 42 
CONTINUE 

WRITE (6 .902) SPLIT .ITER 
FORMAT (/7X.7HSPLIT =»F10.S/ 

7X.36HNUMRER OF ITERATIONS TO FIND SPLIT s»I4//) 


tINl?60 

5.IN1270 

SIN12S0 

SIN1290 

SIN1300 

S.IN1310 

SIN1320 

$IN)330 

SIN1340 

IIN1350 

SIN1360 

SIN1370 

SIN13BQ 

SINI390 

SIN1400 

SIN1410 

SIN1420 

SIN1430 

SIN1440 

$IN1450 

SIN146Q 

IIN1470 

SIN1480 


c 

COMPUTE OMG DIFFERENCE ARRAYS SI... S5 

5IN1490 

c 


SINISOO 


DO 20 J= 2tNN 

<• SINISIO 


JP = J*1 

SIN1520 


JM = J-1 

SIN1530 


SI (J)= PSI ( JP) - PSI (JM> 

$IN1540 


S2< J) = PSI ( JP) - PSI { J) 

ilN1550 


S3( J) = PSI ( J) - PSI ( JM) 

SIN1560 


S4(J) = S2( J)/S3(J)/S1 ( J) 

JIN1570 


S5(J) = S3( J)/S2( J)/S1 (J> 

$IN1580 

20 

CONTINUE 

SIN1590 

C 


SIN1600 

C 

COMPUTE ANGLES THETA (I) 

SIN1610 

C 


S1N1620 


no 410 J=2.NN 

SIN1630 


IF ( (J ,GE. LI). and. (J .LE. L2)) THETA(J*1)= TnOZ 

SIN1640 


IF ( J.LT.Ll ) THETA (J» 1)=THETA( J-1 »1 )-<DN(J)-DN(J-1) )/RCl 

f 1N1G50 


IF ( J.GT.L2) theta (J»1)=THETA( J-1. 1)-(DN(J) -DN( J-1 ) )/RC2 

SIN1660 

410 

CONTINUE 

SIN1670 


DO 2 J=2»NN 

SINI680 


Jv=J-l 

$IN1690 


Da,5*(THETA(J,l)*THFTA(jM.l ) ) 

$IN1700 


X( J)=X(JM) -(nN( J)-DN(JM) )*SIN(0> 

SIN1710 


Y(J) = Y(JM) ♦(riN(J)-DN(jM) )*COS(D) 

$IN1720 

2 

CONTINUE 

SIN1730 


return 

SIN1740 

C . •- : 


IIN1750 

C 

THIS SPLMAX GIVES A TS EXCEEDING TMAX 

5IN1760 

C 

SO TRY LOWERING SPLMAX 1 PERCENT 

$IN1770 

C 


f IN1780 

589 

CONTINUE 

SIN1790 


SPLMAX® SPLMAX * 0.99 

SIN1800 


GO TO 555 

SIN1810 

C 


S1N1820 

C 

NO SOLUTION FOR SPLIT 

SIN1830 

C 


SIN1840 

599 

CONTINUE 

SIN1850 


WRITE(6.9U ) 

SIN1860 

911 

FORMAT (/7X.21HN0 SOLUTION FOR SPLIT) 

SIN1870 


STOP 

SIN1680 


END 

S1N1890 
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1 


c 

c 

c 

c 

c 

c 


c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 


c 

c 

c 

c 


595 


C 

c 

c 

c 


subroutine TMPSFT (N*NN»Rr,fGC»G*DN?*Rl *R^*L1*L2»DNI tONOf 

PSI» R, M* TS» PS* PO* RHO* U» DN* MAXERR) 

Sets an initial ts in middle of nozzle* 

determines corresponding M(D* TS(I)* ON(I)* etc values 
CHECKS difference BETwEEN 0N<N> and DN2 


REAL M(190) 

dimension P0(190) ,THETA(19o*2) *R(190*2) *U(190*2> *151190*2) * 

1 RH0(19U*2) 

dimension OMG(190) *0N(190) ,PS(190*2) 

DOUBLE OMG(190> *DN{i90),PS<190*2> 
dimension PSI (180) 

DOUBLE PSI (190) 

Common /const/ poi*po2*po3* toi*to2*to3* xi* xc, yc 

The starting iteration point is at present determined by the 
location of the nozzle in the duct 
IF nozzle is close to upper wall* the tmax calculation is 
valid only for T5(L1*D* and so nmid should equal LI 
IF NOZZLE in middle OF OUcT NMID SHOULD = (Ll ♦ L2) / 2 
if NOZZLE close TO LOWER WALL NMID SHOULD a L2 
AND CHANGE COMPUTATION OF R(J*NMIO)* AND DN(NMID) 


IF (DNI ,LT. 
IF (DNI ,GE. 
MAXEPR= 0 
Dial ,/(G-l.) 
G2aG/ (G-1 .) 
ROECAYaA. 


DNO) 

DNO) 


NZONEa 0 
NZONE* 1 


Tmax is the temperature which prevents velocity from 

FALLING BELOW .01 (lE GOING NEGATIVE) 

TmAX= (T02/T01 - .0001*( (G-1.0)/(6*2.0) ) ) 

* (P03/P01) «* <(G-U0)/G) * (T01/T02) 

CONTINUE 
TSAa. TMAX 

TSBa TSA - (5.0/TOl) 

0IFFB= 1.0 

Do LOOP WHICH SOLVES FOR TSB SUCH THAT DN2-DN(N)= 0.0 
GUESS NEW TSB USING FIXED TSA AND PREVIOUS TsB VALUE 


STM 

JTM 

STM 

STM 

STM 

STM 

STM 

STM 

STM 

STM 


0 

10 

20 

30 

AO 

so 

60 

70 

80 

90 


STM 100 
STM 110 
STM 120 
STM 130 
STM lAO 
STM 150 
STM 160 
STM 170 
STM 180 
STM 190 
STM 200 
STM 210 
STM 220 
STM 230 
STM 2A0 
STM 250 
STM 260 
STM 270 
STM 280 
STM 290 
STM 300 



IF 

(NZONE 

.EQ. 

0) 

NMI0= L2 





STM 

310 


IF 

(NZONE 

.EQ. 

0) 

R(NMIO,I)a 

R1 

« 

EXP( 

-ROECAY*r1*(DNI»1.0) ) 

STM 

320 

1 






R2‘ 

« 

EXP( 

ROECAY*P2*ONO) 

STM 

330 


IF 

(NZONE 

.Eo. 

0) 

DN(NMIO)= 

ONI 


1.0. 


STM 

3A0 


IF 

(NZONE 

• EO. 

I ) 

NMID= Ll 





STM 

350 


IF 

(NZONE 

.EO. 

1) 

R(NMlO*l)a 

R1 

« 

EXP( 

-RDECAY*R1*DNI) 

STM 

360 

I 






♦ R2 


EXP( 

RD£CAY*r2*(ONO*1.0) ) 

STM 

370 


If 

(NZONF. 

.EO. 

1) 

DN(NMlD)a 

DNI 




STm 

380 


STM 390 
STM AOO 
STM AlO 
STM A20 
STM A30 
STM AAO 
STM AA5 
STM A50 
STM A60 
STM A70 
STM A80 
STM A90 
STM 500 
STM 510 


DO 

AO J=l*100 


STM 

520 

IF 

(J .£0. 1) 

TS(NMID*l) = TSA 

STm 

530 

IF 

(J .GE. 3) 

TsBa (DIFF8«TSA-DIFFA*TSB) / (DIFFB-DIFFA) 

STm 

5A0 

IF 

(J .6E. 2) 

TS(NMZD,l>a TSS 

STM 

550 

IF 

(TS(NMI0*1) 

•GT. TMAX) GO TO 599 

STM 

560 
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C SOLVE FOR TEMPfRaTURESvMACH NUMBFRS*ETC, AT NODE POINTS 

C WORKING FROM NOZZLE MIDDLE (NMID) OUT TO BOTH WALLS 

C 

c nozzle centre line values 

c 

M(NMI0)= SORT-(2.0 * G1 (1.0/TS(NMID*1> - loO) > 

PS<NMID»1)= TS(NMID*I> G2 
ftHO(NMIDfl)= PS(NHID*1) / TS(NMlD9l) 

U(NMIO*I)= MCNMID) « SQRT<G » TS^nMIOjD) 

P0(NMID)= PSINMID,!) « <1,0 ♦ 0o5 / Gl« M{NMID)«*2) *»62 
C 

C REGION EROM NOZZlE MIDDLE TO OUTER WALL 

C SET M* TS* PS 

C <L2 IS JUMP FROM NOZZLE TO SECONDARY STREAM) 

C 

PK=1. 

OK=l, 

1 

NPRR=NMID+1 
DO 42 I=NPRRfN 
IM=I-1 
C 

C JUMP POINT TO SECONDARY STREAM 

C 

IF (I .E0.L2+1) PK= P02/P01 

IF (I .EO.L2+1) QK=P02/P01»(T01/T02)»»G2 

If (I .EO.L2+I) RK=T02/T01 

SX=RK»R(IM»1)*(PSI (I)-PSI (IM) ) 

M(I>= M(IM) - SX/{G*«.5 * QK e TS ( IM» 1 ) »» <G2->0 o5) ) 

Bl= ((PK/PS( IH*1)> «»(lo/G2)-lo)»2.«Gl 
IF( (I.E0.L2^1 ) .and, (B l.LEoOoO) ) TMAX=TMAX*0 .995 
IF( (I.E0.L2+1) .AND. (BKLE.O.O) ) GO TO 595 
IF (1 .EG. L2*l) M(I>= SQRT(bI) 

C 

C set TS* PS* RHO* U* P0» Dn* R ACROSS STREAM FROM M(I> 

C 

TS(I*1)=RK/(1,*0.5/61«(M(I) )»»2o) 

PS { I * I ) =OK" < TS ( I ♦ I ) ) «<>G2 
RHO { I * 1 ) sPS (I * 1 ) /TS ( I * 1) 

U ( I , 1 ) =M ( I ) »SORT (G«TS (1*1)) 

PO tI)=PSa*l)»a.*0.5/Gl*(M(l) )»»2,)*»G2 
Z= 2.0 / {RH0(lM*l)«U{lM*l>+RH0(l9l>*Ua9l) ) 

On(I)=: 0N(1M) + Z » (PSI(I) - PSI(IM)) 

Dwl= DN(I) 

0W2= 0N2 - DN<I) 

IF tDW2 .LT. 0,0) 0W2= 0.0 

42 R ( I * 1 ) =R1 * ( 1 . iDWl/DN2) «EXP (-RDECAY«DW1*A8S (R1 ) ) «R2«DW1/DN2» 

I EXP(-RDECAY«DW2*ABS(R2) ) 

C 

C search FOR TSA WHERE 0N2-0N(N) .EQo 0,0 

C 

IF (J .EG, 1) DiFFAs DN2-DN(N) 

IF {J .6E, 2) DTFFBs DN?-DN(N) 

C 

C JF SUFFICIENTLY SMALL INTERVAL? THEN EXIT 

C 

IF lABS(TSB-TSA) ,LE. 0,0000001) GO TO 41 
IF (ABS(0IFF8) .LE. O.OOOl) GO TO 41 

40 continue 

WRITE(6*900) DN2*0N(N) *TSa*TSB 

900 FORMAT (/7X*22HN0 CONVERGENCE DN2 =*f10.5,5X*7hDN<N) =*f10<,5*5X 
1 5HTSA =*F10.5*5X*5HTSB -9F10,5/> 

41 CONTINUE 


STM 570 
STM 580 
$TM 590 
ITM 600 
STM 610 
STM 620 
STM 630 
STM 640 
STM 650 
STM 660 
STM 670 
STH 680 
STM 690 
STM 700 
STM 710 
STM 720 
STM 730 
STM 740 
STM 750 
STM 760 
STH 770 
STM 780 
STM 790 
$TM 800 
STM 810 
STM 820 
STM 830 
STH 840 
STM 850 
$TM 860 
$Tm 870 
STM 872 
STM 873 
STM 880 
STM 890 
STM 900 
STH 910 
STM 920 
STH 930 
STM 940 
STM 950 
STM 960 
STM 970 
STM 980 
STM 990 
STMIOOO 
STMlOlO 
STM1020 
STM1030 
STM1040 
STM1050 
$TH1060 
$TH1070 
STM1080 
STM 1090 
STM 1100 
STMlllO 
STM1120 
STH 1130 
STM I 140 
STM 1150 
STM 1160 
STM1170 
STMllSO 
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ono ooo ooo 


kegion from nozzle middle to inner duct wall 


PK=1, 

0K=1. 

RK=1. 

NPRR=NMI0-1 
DO 46 KK5=1*NPRR 
I=NMID-KK 
IM=I+1 


JUMP point to tertiary stream 

IFU ,EQ. Ll-1) PK=P03/P01 
IF(I *EO. Ll-1) OK = P03/P01*(T01/T03)«<»G2 
JF(I .EO. Ll-I) RK 5 T03/T01 
SX=RK*RaM«l)*{PSl{I)-PSI (IH) ) 

Mai- MUM) - SX/tG**.5 OK * TS( IM, 1 ) *» { G2*0.5> > 

61= ((PK/PS( IM*1>) *»(1,/G2)-1.)*2.*G1 

IF<<1,£Q,L1*1? (BI.LE.0,0) ) TMAX=TMAX*0#995 

IFUI, EQ, Ll-1) .AND, tBl.LE. 0.0)) GO TO S95 
IF (I .£0. Ll-1) M(I)=SQrT(B1) 

SET TS* PSf RHO* U, PO* ON» R ACROSS STREAM FROM MH) 

TS(I»l)=RK/a. + O.S/Gl*(H{l) )**2.) 

PS(itl)=0K«tTSa,l))**G2 
RHO (1*1) -PS (I » 1 ) /TS ( I » 1 ) 

U ( 1 , 1 ) =M ( I ) <»SQRT (6»TS ( I * 1 ) ) 

PO n ) sPS ( I ♦ I ) « a , ♦ 0 .5/Gl« (M a > ) ««2 , > *«G2 
Z= 2.0 / (RHOtIM*l)»U(lM,l)+RHO(I,l)*U(I>l>) 

DN(I)= ON(XM) * Z * (PSI(I) - PSKIM)) 

DW1= DN(I> 

IF (DWl .LT. 0.0) DW1= 0.0 

DW2= DN2 - ONtI) ' 

46 R a » 1 ) =Rl* U ."i-DWl/ONa) *EXP (-ROeCAY*Owi« ABS (R1 ) ) ♦R2«0W1/DN2* 

1 EXP(^ftO£CAY*DWZ«ABS{R2) ) 

RETURN 

599 CONTINUE 
MAXERRs I 
RETURN 
END 


STM1190 

$TM1200 

$TM1210 

STM1220 

5TM1230 

$TK1240 

STM1250 

STM1260 

STM1270 

5TMJ280 

STM1290 

STM1300 

STM1310 

5TM1320 

STM1330 

STM1340 

$TM1350 

STM1360 

STM1370 

$TM1373 

$TM1374 

STM1380 

5TM1390 

STM1400 

STM1410 

STM1420 

STM1^30 

$Tm1440 

$TM1450 

STM1460 

STM1470 

STM1460 

$Tm1490 

STM1500 

STMISIO 

STM1520 

STM1530 

STM1S40 

STM1550 

$Tm1560 

5TM1570 

sTMisao 
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20 


«;UBR0UTIN£ EDOY (N»NN»NX*UtPS*CC»OS»E»El»KHOtVlS«R* Y* 
1 S4f S5* OELl,OEL2tDEL3,OEL4,COREl*CORe2»CORE3) 

«•* *•»«■ *■»«■»<>■*»■»*«■*»»*«•«•**#*********»»** * 

SUBROUTINE EOnY 

CALCULATES EDDY VISCOSITY VALUES - ARRAY E(N) 

REQUIRED PARA?^ETEPS: 

ARRAYS 0MG» RHO> ViS, Yt PS> DS* Sl.*.S5 
VARIABLES Nt NN* N1...N6* COREl* COREPt C0RE3 
METHOD 

CALCULATES REFERENCE VARIABLES (DP* DELI* DELO* 

UREFI* UREFO* CC) 

FOR EACH J FROM 1 TO M 

DETERMINES THE REGION BY CHECKING N1-**N5 
CALCULATES EDDY VISCOSITY: E(J) 


REAL LZ*L7R*LM 

DIMENSION SI U90) *S? ( 190) *S3 (190) *S4 (190) *S5 ( 190) *E(190) *E1 (190) * 

1 OS (190)* VIS (190) *THETA(l90*2)»R(l90f2)*U(190*2)* 

2 TSU90*2> *RH0(19Ct2) *NX(6) 

DIMENSION OMG(190)» Y ( 190) *PS ( 190 »2) 

DOUBLE OMG(190)> Y(190.> *PS(190*2) 

CALCULATE reference VARIABLES 

np=(PS(i*2)-PS(i.l) )/os(i) 

DELI=Y(2)-Y(1) 

DELO=Y(N)-Y(NN) 

UREFI=SORT( (U( 2*2) /DELl-OELI*DP/ (2.*VIS (2) ) ) / RHO( 2*2) * 

1 Vise 2)) 

nP-(PS(N»2)-PS(N*l) )/DS(N) 

UREFO=SORT( (U(NN*2)/DEL0-DEL0*0P/(2.*VIS(NN)> )/RH0(NN*2)* 

1 VIS(NN)) 

IF (COREl .EO. 0.0) CC- 0.08 
SCORE = CORE! ♦ C0RE2 ♦ C0RE3 

DO LOOP WHICH COMPUTES EDDY VALUES E(J) 

YS: DIMENSIONLESS OlST TO WALL 

LZR: curvature EFFECT 

RN: RICHARDSON NUMBER 

DIFFERING CALCULATIONS FOR EACH REGION (N1...N6) 

no 10 J-2*NN 

.JM-J-l 

JP=J*1 

IF (J. .LT. NX(l) ) GO TO 20 
IF (J .GT. NX (61 ) 60 TO 21 

IF (J .LE. NX(2) ) 60 TO 30 

IF (J .GE. NX (5) ) GO TO 30 

IF (J .LT. NX(3) ) GO TO 29 

IF (J .GT. NX(4) ) GO TO AO 

GO TO 30 

INNER WALL 

DRN=Y(J)-Y(1) 


SED 0 
SEO 10 
JED 20 
SED 30 
$ED AO 
SED 50 
$ED 60 
SED 70 
SED 80 
SED 90 
SED 100 
SED 110 
SED 120 
SED 130 
SED lAO 
SED 150 
SED 160 
SED 170 
SED 180 
SED 190 
SED 200 
SED 210 
SED 220 
SED 230 
SED 2AQ 
SED 250 
SED 260 
SED 270 
SED 280 
SED 290 
SED 300 
SED 310 
SED 320 
SED 330 
SED 3A0 
SED 350 
SED 360 
SED 370 
SED 380 
SED 390 
SED AOO 
SED AlO 
SED A20 
SED A30 
SED AAO 
SED ASO 
SED A60 
SED A70 
SED ASO 
SED A90 
SED 500 
SED 510 
SED 520 
SED 530 
SED 5A0 
SED 550 
SED 560 
-SED 570 
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YS=bHN*UREri«RHO( J»1)/VIS(J> 

LZ=.41»(1.-EXP(-YS/26,) )*nRN 
AZ1=0*C9*DEL1 

nD2= 1.0 - 2.0 » DRN / Y{N> 

TF (SCORE. EO. 0.0) AZ1== ( . 1A-.08*DD2*«2-.06^»DDZ**^) «Y (N) » ,5 
LZ= LZ * 1.2 
IFtLZ .67. AZl) LZ= AZl 
^ 00 TO 25 

C 

C CORE REGIONS 

C 

30 E(J)=0. 

r-0 TO 10 
C 

C DEL2 REGION 

C 

29 NX3=NX(3) 

NX2=NX(2) 

LZ=CC*DEL2* ( 1 •♦.6*U (NX2»2) /U (NX3»2) ) 

GO TO 25 
C 

C 0EL3 REGION 

C 

40 NX5=NX(5) 

NX4=NX(4) 

LZ=QC<^DEL3*(l,*.6»U(NX5*2)/U{NX4t2)) 

GO TO 25 
C 

C OUTER WALL 

C 

21 CONTINUE 

nRN= Y(N)-Y(J) 
ys=drn«urefo«rho ( j* 1 ) /vis ( J) 

LZ=C.41*(l.-EXP(-YS/26.) )«DRN 
aZ4=0.09«DEL4 
DDZ=1 . 0-2. 0*DRN/Y (N) 

IF (SCORE. EQ. 0.0) AZ4=(0. 14-0.08»D02**2-0.06*ODZ**4)*Y(N)*0.5 
IF(LZ .GT. A24) LZ=AZ4 
LZ= L2 * 1.2 
25 CONTINUE 

PN= 2.0 * (R(J*1> / RHO(J*l)> 

PN=RN/(S5(J)»(U(JP*1)-U(J,1) >*S4(J>«(U(J»1)-U{JM,1^>) 

IF (RN .LE. 0.0> LZR= LZ*(2.0 - EXP(3.0*RN)> 

IF (RN .GT. 0.0) LZR= LZ*EXP (-3«0*RN> 
el=Y(J)-Y(JH) 

EJ=Y(JP)-Y(J) 

EK=Y(JP)-Y(JM) 

DUY=AB3(EI*(U{JP*2)-U(J*2) )/(EK*EJ)+EJ*(U(J*2>-U(jM,2))/(EK*ED) 
E(J)=DUY*LZR«LZP ' ' ' 

10 CONTINUE 

E(l)= 0.0 
F(N)= 0.0 
RETURN 

eno 

SUBROUTINE WALLS (XW,YW, XX, YY*T»CURtM*MQ) 

C subroutine walls smooths THF boundary data USING A LEAST 

C SOUARES PROCEDURE. IT ALSO INTERPOLATES THE SMOOTHED 

' c Data to get the curvature and slope at any point 

Q «««««««»««« 

dimension XW(99,2) *YW(99*2) *YP(99*2) *CURV(99,2) *NOl (10#2) * 

1 YB(99,2) 


$ED 5B0 
tED 590 
SEO 600 
SED 610 
SEO 620 
SEO 630 
SEO 640 
SED 650 
SED 660 
SED 670 
SEO 680 
SED 690 
SEO 700 
SED 710 
SEP 720 
SED 730 
SED 740 
SEO 750 
SED 760 
SED 770 
SED 780 
SED 790 
SED 800 
SED 810 
SED 820 
SEO 830 
SED 840 
SED 850 
SED 860 
SEO 870 
SED 880 
SED 890 
SED 900 
SED 910 
SED 920 
SED 930 
SED 940 
SED 950 
SED 960 
SED 970 
SED 980 
SED 990 
SEDIOOO 
SEOlOlO 
SED1020 
SED1030 
SED1040 
SED1050 
SED1060 
SED1070 
SED1080 
SED1090 
SEDllOO 
SEDino 
SWA 0 
SWA 10 
SWA 20 
SWA 30 
SWA 40 
SWA 50 
SWA 60 
SWA 70 
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COMMON/KUHV/MST *MSS*NOI fCllRV 

»WA 

80 


TF(MO.EO.l) MsMST 

SWA 

90 


IF(M0,EQ.2) M=MSS 

SWA 

100 


IF<XX.LT,999.) GO TO 60 

SWA 

110 


ri=N+M-i 

SWA 

120 


MX=N 

SWA 

130 


MP=M-1 

SWA 

140 


00 4 Jsl.MP 

SWA 

150 


JK=M-J 

SWA 

160 


JKlsJK*! 

SWA 

170 


Nt>=NDl (JKl »MQ)-NDI ( JKfMQ)*l 

SWA 

180 


00 3 IsltNP 

SWA 

190 


K=NX+1-I 

SWA 

200 


l=k-jk 

SWA 

210 


XW(KfMO>sXW(L*HO) 

SWA 

220 

3 

YW(KfMO)=YW(L*MQ) 

SWA 

230 


NOI < JKl *MO)=NX 

SWA 

240 

A 

NX=NX-NP 

SWA 

250 


DO 5 I=lfN 

SWA 

260 

5 

YBn*MO)=YW<I*MQ) 

SWA 

270 


KM = 1 

SWA 

280 


no 198 MPsl*M 

SWA 

290 


N0A=NDI <MPtHO> 

SWA 

300 


kl=nda-km*i 

SWA 

310 


TF(KL.LT,3) GO TO 700 

SWA 

320 


1F(KL.EQ.3> go to 30 

SWA 

330 


TF(KL.EQ.4> go to 40 

SWA 

340 


L=5 

SWA 

350 

16 

no 20 KS=KM*NDft 

SWA 

360 


IF<L.E0.4) GO TO 17 

SWA 

370 


CALL L50(KS»K>ltN0AtL»Al«A2«A3tXW(l«MQ) tYBOvNOl?- 

SWA 

380 

17 

YH <KS»MO>=AI*A2»XW(KS»MQ)^a3*XW(KS*MQ)**2 

SWA 

390 


YP (K5»M0) =A2*2,*A3*XW (KS*MQ> 

SWA 

400 


CURV (KS.MO) =-2.*A3/U •♦YP(KS»MQ)**2)**1,5 

SWA 

410 


TF (APS (CURV (KS*MO) ) -LT..001) CURV < KS* MO) *0, 

SWA 

420 

20 

CONTINUE 

SWA 

430 


GO TO 198 

SWA 

440 

30 

C1=XW (KM,MO)-XW (KM*1 tMO) 

SWA 

450 


C?=XW (KM*MQ) -XW (NO A, MO) 

SWA 

460 


C3=XW(KM*1 ,M0)-XW(NDA,M0) 

SWA 

470 


YPP=2,* (YW (KM*MO) /C1/C2-YW (KM o 1 • MQ ) /Cl /C3« YW (N0 A« mO) /C2/C3) 

SWA 

480 


no 31 K=KM*NnA 

SWA 

490 


YP (KfMO)= (2,*XW (K*MO)-XW<NnA*MQ)-XW(KM+l *M0) ) *YW (KH«MQ) /C1/C2 

SWA 

500 


1 -<2.*XW(K«M0)-XW(KM*M0)-XW (mDA*MO) ) *YW (KM« 1 «MQ) /C1/C3 

SWA 

510 


2 ♦ (2.*XW(K*M0)-XW<KM,M0)-XW(kM*1 *MQ> )*YW<N0A*MQ)/C2/C3 

SWA 

520 


CL»PV (K*MQ) =-2.*YPP/( 1. + YP(K*M0)**2)**1.5 

SWA 

530 


IF(ABS(CURV(K.M0) ).LT,0.001> CURV (K.MO) *0. 

SWA 

540 

31 

CONTINUE 

SWA 

550 


OO TO *198 

SWA 

560 

40 

L=4 

SWA 

570 


CALL LSQ(KS9KMfN0A»LtAl«A?«A3*XW<l«MQ) tYBUtMO) ) 

SWA 

580 


GO TO 16 

SWA 

590 

198 

KM=NDA*1 

SWA 

600 


return 

SWA 

610 

60 

NDA=1 

SWA 

620 


no 70 I=1»M 

SWA 

630 


J=N0I (I*MO) 

SWA 

640 


IF((XX,GE.XW(NDA«H0)).AN0.(XX.LE.XW(J»HQ>)) go TO 71 

SWA 

650 


IF(I.EO.M) go to 71 

SWA 

660 

70 

N0A=J*1 

SWA 

670 

71 

JJ=J-1 

SWA 

680 


00 72 K*NDA*JJ 


89 


72 

73 


700 

900 


4 

5 


C 

c 




SWA 

690 

lF<(XX.GE.XW(K,MO)).ANr).(XX,LE*XW(K*l»MO))> GO TO 

73 

SWA 

700 

CONTINUE 


SWA 

710 

IF(K.EQ.NDA) K=NnA+l 


SWA 

720 

C1=XW (K>1 *MQ)-XW (K*MQ) 


SWA 

730 

C?=XW<K*l*MO)-XW(K-l *MO) 


SWA 

740 

C3=XW(K*MQ)-XW(K-1 *MQ) 


SWA 

750 

SlsXX-XW(K+l»MQ) 


SWA 

760 

S?=XX-XW(KtMO> 


SWA 

770 

S3=XX-XW(K-1»M0) 


SWA 

780 

YY=S3*S2«YW<K+1 *MO)/Cl/C2-Sl*S3*YW{K»MO)/Cl/C3*Sl*S2*YW(K-l*MO>/ 

SWA 

790 

1 C2/C3 


SWA 

800 

CUPsS2*S3*CURV(K+lfMO)/Cl/C2-Sl*S3*CURY(KfMO)/Cl/C3 


SWA 

810 

1 ♦Sl*S2«CURV(K-l»MQ)/C2/C3 


SWA 

820 

T = S3*S2*YP(K*l,M0)/Cl/C2->Sl*S3*YP<KtMQ)/Cl/C3*Sl*S2*YP(K-l*M0> 

SWA 

830 

1 /C2/C3 


SWA 

840 

TsATAN (T) 


SWA 

850 

RFTURN 


SWA 

860 

WRITE (6»900) 


SWA 

870 

Format (//4X*34H less than three points in segment) 


SWA 

880 

STOP 


SWA 

890 

END 


SWA 

900 

SUBROUTINE USD (KS *KM*NDA»L ♦ A1 ♦ A2t A3*Ct8) 


SLS 

0 

DIMENSION B(99> *C(99) *S(2*4> 


SLS 

10 

K = KS 


SLS 

20 

IF ( (KS.EO.NDA) .OR, (KS.EQ.NDA-1 ) ) KsNDA-2 


SLS 

30 

IF( (KS.EO.KM) .OR. (KS.EQ.KM*!) ) K«KM*2 


SLS 

40 

IF(L.E0.4) KrKM*? 


SLS 

50 

no 5 I=l»2 


SLS 

60 

r«0 S J=1.4 


SLS 

70 

M=I-1 


SLS 

80 

S(I*J)=0. 


SLS 

90 

no 4 M=1,L 


SLS 

100 

KK=K*M-3 


SLS 

110 

TF(C(KK) .EQ.O.) C(KK)=l.E-06 


SLS 

120 

IF(B(KK) .EQ.O.) B(KK)=l.E-06 


SLS 

130 

S(I»J)= S(I*J) *B(KK)**N*C(KK)**< J-N) 


SLS 

140 

CONTINUE 


SLS 

150 

A=L 


SLS 

160 

D=A»(S(1.2)*S(1 »4)-S(l t3)**2)-S(l.l)*(S(l»!)*S(l,4) 

-S(1,2)*S(1,3) 

SLS 

170 

1 ) *S( 1»2>*(S(1.1)*S(1.3)-S<1»2)**2) 


SLS 

180 

A1-(S(2»1)*(S(1*2)*S(1,4)-S<1 ,3)**?)-S(2*2)*(S(l,l)*S(1.4)- 

SLS 

190 

1 S(1.2)*S(1»3) )*S(2»3)*(S( l*l)*S(l,3)-S(l»2)**2) )/D 


SLS 

200 

A2=(S(2tl)*(S(l*2>*S(l»3)-S(l,l)*S(l»4) ) *S ( 2 , 2) * ( A*S ( 1 .4) - 

SLS 

210 

1 S(1.2)**2)*S(2»3)*(S(ltl)*S(l»2)-A*S(l»3)) ) /O 


SLS 

220 

A3=(S(2*l)*(S(l»l)*S(l,3)-Sn,2)**2)-S(2*2)*(A*S(l* 

3>- 

SLS 

230 

1 S(1.1)*S(1»?) )*S(2*3)*(A*S(1,2)'-S(l*l>**2) )/0 


SLS 

240 

RETURN 


SLS 

250 

END 


SLS 

260 

SUBROUTINE SOLY (OStN.NNfPPtPRTfGI- 


SSO 

0 

REAL M(190) 


5S0 

10 

DIMENSION SK190) .52(190) » S3 (190) »S4< 190) *S5 (190) ♦£( 190) »E! (190) f 

SSO 

20 

1 DS(190) »VIS(190) ,X(190) »Y(190)*P0(190) »H(190) ♦ 

SSO 

30 

2 THETA (190,2) ,R ( l90 .2) ,U ( 190 ,2) ,TS (190 ,2) »RHO( 190,2) , 

SSO 

40 

3 A(190) ,R(190) .C(190) ,0(190) 


SSO 

50 

DIMENSION OMG(I90) ,ON(190) ,PS(190,2) 


SSO 

60 

DOUBLE OMG(190) ,DN ( 1 90 ) , PS ( 1 90 , 2) 


SSO 

70 

DIMENSION PSI (IBO) 


SSO 

80 

DOUBLE PSI(190) 


SSO 

90 

COMMON/ INCD/ U.PS,TS,M,PO,TH£TA,R,X,Y,ON 


SSO 

100 

COMMON/SOLV/ SI ,S2,S3,E,El ,H,S4,S5 


SSO 

110 

COMMON/INEOY/PSI,RhO,VIS 


SSO 

120 

CL=(G-1.)/G 


SSO 

130 


CO « 0. 


90 



oooo ooooooooo 


IF(PR .EO, .7) CO = 1* 

Do 10 J=2*NN 

JP=J*1 

JM= J-1 

Y1 =E1 tJP)*<VIS< JP)/PR ♦ RhO(JP*1)*F<JP)/PRT> 

Y2sEl(JM)*( VIS(JM)/PR ♦ RhO < JM* 1 ) *E < JM) /PRT) 

Y? »E1(J) *(VIS(J )/PR ♦ RH0(J *1)*E(J )/PRT) 
yt* » (Y1*Y3)/S2<J) 

YS s<Y3^Y2)/S3(J) 

Visa Y4 / SKJ) 

V7s Y5 / SKJ) 

Y8= Y6 • U(J*l) 

Y9= Y7 * U(J,l) 

A(JM) * U(J*l)/OS<J) ♦ Ya ♦ Y9 
B(JM) s -Y8 
C<JM) * -Y9 

IF (PR. NE. 0.7) D(JM)=U(J*1)**2/DS(J)-(PS(J*2)-PS(J*1> >/DS(J)*(l./ 

1 (RHO(J*n)> 

TF(PR.E0.0.7) n(JM)*TS(J*l )*U(J.l)/DS(J)*(CL*U(J,l)/(RHO(J*l) ) > 

1 *(PS(J*2)-PS(Jf 1) )/DS( J)«RH0(J«1>*CL*(VIS(J) «RHO( Jf 1)*E( J) >*(U 

2 ( J»1 )**2*(S5(J)*(U(JPt2)-U(J*2) )♦S4(J)*(U(J♦2)-U( JM*2) ) > )**2 

10 CONTINUE 

A(l) sA(l) * C0»C(1> 
a(N-2) = A(N-2) ♦ C0*B(N-2) 

CALL CALC(A«B*CfD«H*N) 

413 RETURN 
EMD 

Subroutine CALC(A«BfC«o*Hf j> 

DIMENSION A(190) 9B(190)«C(190) »0<190> •W(190) «(><190> fH(190) •QU90) 
K2 s J-2 
Ml= J-2 
W(l) = A(l) 

G(l) s 0(1) /W (I) 

Do IKS ?,N2 
K1 s K-1 

0(K1) s R(Kl)/W(Kl) 

W(K) s A(K)-C(K)*Q(K1) 

1 r,(K) s ( D(K)-C(K>*G(K1) )/W(K) 

H(N2) s G(N2) 

N3 =J-3 

no 2 K * 1*N3 

KK = N2-K 

2 H(KK) s G(KK) -Q(KK)*H(KK*1) 

return 

End 

subroutine DSMOVE(XW«Ytft XX fYYvOSXt theta tRR«NPAIR,J) 


SUBROUTINE DSMOVE 

FINOS NEW X VALUE AT A GlvEN DISTANCE DSX ALONG WALL SURFACE 

method: 

GUESSES VALUES OF DX UNTIL (DSX - DS COMPUTED) s Q.O 


DIMENSION XW(99*?> .YW(99*2) 

APPROXIMATE POINT X(N)* Y(N) 

AT distance DS(N) ALONG UPPER WALL 
(TWO METHODS OF APPROX IMATING* DEPENDING ON Rl) 


SSO 

140 

$so 

150 

SSO 

160 

ISO 

170 

ISO 

180 

SSO 

190 

ISO 

200 

ISO 

210 

SSO 

220 

ISO 

230 

SSO 

240 

ISO 

250 

ISO 

260 

ISO 

270 

SSO 

280 

SSO 

290 

ISO 

300 

ISO 

310 

SSO 

320 

ISO 

330 

SSO 

340 

SSO 

350 

SSO 

360 

ISO 

370 

SSO 

380 

SSO 

390 

SSO 

400 

ISO 

410 

ICA 

0 

SCA 

10 

SCA 

20 

SCA 

30 

ICA 

40 

SCA 

50 

SCA 

60 

SCA 

70 

SCA 

80 

SCA 

90 

SCA 

100 

SCA 

110 

SCA 

120 

SCA 

130 

SCA 

140 

ICA 

150 

SCA 

160 

SCA 

170 

IDS 

0 

IDS 

10 

SOS 

20 

IDS 

30 

SOS 

40 

IDS 

SO 

IDS 

60 

IDS 

70 

IDS 

80 

IDS 

90 

SOS 

100 

IDS 

110 

SOS 

120 

IDS 

130 

SDS 

140 


91 


ooooooooooooooooooo oooo 


510 

900 

512 


500 


OXST/VRs DSX * ARS{COS(THeTA) ) 

XSTAOs nXSTAR ♦ XX 

CALL WALLS<XW,YW, XSTAR* YSTAR» TSTAR* RSTARfNPAjR. J> 

DDSTARs DSX - CURVDS(DXSTAR, RR»RSTAR» THFTA»TSTAR*YY* YSTAR) 
nXBs DSX * ABS(CnS( (THETA ♦ TSTAR)*0,5)) 

Xfl- DXB ♦ XX 

CALL WALLStXW,YW,XB* YB* TB* RB» NPAIR» J) 

DOIFFB= DSX - CIJRVDS(DX8, RRt RB, THETA* TB.YY* YB) 
no 610 1=1*50 

IF (ABS(DDIFFB-DDSTAR) ,LE. 0.1 «*8) GO TO 512 
nXB=(DDIFF8*DX5TAR - 0DSTAR*0XB) / (DDIFFB-DDSTAR) 

IF (DXB .LT, 0.0) DXB= DXSTAR * 0.1 
XB= DXB ♦ XX 


CALL WALLS(XW* YW* XB* YB* TB* 
ODlFFBs DSX - CllRVDS(DXB. RR* 

CHECK FOR termination 


RB* 

RB* 


NPAIR* 

THETA* 


J) 

TB*YY* YB) 


60 TO 512 
60 TO 512 


IF (ABS(XB“XSTAR> .LE. 0.1 **6> 

IF (ABS(DDIFFB) .LE. 0.1 ** 7) 

CONTINUE 

WRlTE(6*q00) XB. DDIFFB 

format (/7X.21HN0 CONVERGENCE DSM0VE*SX*4HXB 
1 Flo. 5/) 

CONTINUE 
XX* XB 
RETURN 
Emd 

real function CURVDS(DX*Ra*RB*TA*TB*YA*YB) 

IF (ABS(TA-TB) .LE. 0.0001) GO TO 500 
DR= ABS(RA*RB) 
if (DR .LE. 0.02 ) GO TO 500 

CURVDS* (2.0/DR) * ABS(TB-TA) 

return 

CURVDS= SQRT(DX ** 2 + ( ABS ( YB-YA) **2) ) 

RETURN 

Emd 

SUBROUTINE LOOK (N . NN . K3 * NCORE » DEL 1 * DEL2 * DEL3* DEL4. 

1 COREl .C0RE2. CORES. NX. PO. DN. U* DlPi 

««»•«»•«««***««««»««*»»*»«»«■«««««««»««««««* 

SUBROUTINE LOOK 

CALCULATES 

EDGES OF POTENTIAL FLOWS - Nl... Nti 
WIDTH OF FLOWS - DELI, DEL2. DEL3. DEL4 
width of cores - CORE)* C0RE2, C0RE3 
REQUIRED PARAMETERS: 

ARRAYS DN. PO 

VARIABLES N. NN. NST£P* K3. POl* P02. P03 
METHOD 

CALCULATES TOLERANCES TOLl. TOL2, T0L3. TOLA 
CALCULATES ARRAY DIP • RATE OF CHANGE OF pO 

FOR EACH J FROM 2 TO NN 
SCAN ARRAY DIP(J) 

FOR REGION WHERE DIP GT TOLERANCE 
SET FLOW REGION EDGE NX 
SET FLOW REGION WIDTH DELX 
FOR REGION WHERE OIP(J) LT TOLERANCE 
SET FLOW REGION EDGE NX 
COMPUTE CORE REGION AREA COREX 
♦*♦•****•****#»* **********«•*****■»**•»***#*#** ***♦«*•** 


NMAX) 


IDS 150 
SDS 160 
$DS 170 
$05 180 
SDS 190 
SDS 200 
SDS 210 
SDS 220 
SDS 230 
SDS 240 
SOS 250 
SDS 260 
SDS 270 


'F10.A*SX.8HDDIFF8 ** 


SDS 

280 

SDS 

290 

SDS 

300 

SDS 

310 

SOS 

320 

SDS 

330 

SDS 

340 

SOS 

350 

SDS 

360 

SDS 

370 

SDS 

380 

SDS 

390 

SDS 

400 

SDS 

410 

SDS 

420 

SDS 

430 

SCU 

0 

SCU 

10 

SCU 

20 

SCU 

30 

SCU 

40 

SCU 

SO 

SCU 

60 

SCU 

70 

SCU 

80 

SLO 

0 

SLO 

10 

SLO 

20 

SLO 

30 

SLO 

40 

SLO 

50 

SLO 

60 

SLO 

70 

SLO 

BO 

SLO 

90 

SLO 

100 

SLO 

no 

SLO 

120 

SLO 

130 

SLO 

140 

SLO 

150 

SLO 

160 

SLO 

170 

SLO 

180 

SLO 

190 

SLO 

200 

SLO 

210 

SLO 

220 

SLO 

230 

SLO 

240 


92 



c 



SLO 

250 


dimension X(190> tYdOO) tUM90«2) «DIP(190) »NX(6) •P0<190) ;H(190) 

SLO 

260 


DIMENSION DN<190) 


SLO 

270 

c 

DOUBLE DN(190) 


SLO 

280 


COMMON/CONST/ P01tP02»P03* T01tT02*T03* XI* XC* YC 

SLO 

290 

c 



SLO 

300 

c 

SET TOLERANCES 


SLO 

310 

c 



SLO 

320 


TOL1=0*01*P03/OM<N) / POl 


SLO 

330 


TOL2=0.01*(P01-P03)/ON(N) / POl 


SLO 

340 


T0L3= 0.01*(P01-P02)/DN(N) / POl 


SLO 

350 


T0L4* 0-01*P02/DN(N) / POl 


SLO 

360 


1F(K3,EO.O)CORE3=0* 


SLO 

370 

c 



SLO 

380 

c 

This loop finds a value for nmax 


SLO 

390 

c 



SLO 

400 


DO 11 J*2*N 


SLO 

410 


OlP(J)s(PO(J)-PO(J-l))/(ON(J>-DN(J- 

l>> 

SLO 

420 


OlP(J)= ABS(OIP(J)) 


SLO 

430 

11 

Continue 


SLO 

440 


NMAXsl 


SLO 

450 


DO 12 J=?»N 


SLO 

460 


IF (UIJ*2> .GT. U<NMAX*2) ) NMAX* 

J 

SLO 

470 

12 

CONTINUE 


SLO 

480 

C 



SLO 

490 

c 

THIS LOOP COMPUTES EDGE VALUES FOR THE PRIMARY, .SECONDARY 

SLO 

500 

c 

AND tertiary potential 

FLOWS 

SLO 

510 

c 

tertiary FLOW REGION 


SLO 

520 

c 



SLO 

530 


DO 50 I=2*N 


SLO 

540 


IF (OIP(I) ,LT. TOLU GO TO 52 


SLO 

550 

50 

CONTINUE 


SLO 

560 

52 

CONTINUE 


SLO 

570 


NX(1>= I-l 


SLO 

580 

C 



SLO 

590 

C 

SECONDARY FLOW REGION 


SLO 

600 

C 



SLO 

610 


DO 60 M=2.N 


SLO 

620 


T= N ♦ (2-M) 


SLO 

630 


IF (DIP(I) ,LT. TOLA) GO TO 62 


SLO 

640 

60 

CONTINUE 


SLO 

650 

62 

CONTINUE 


SLO 

660 


NX(6>= I 


SLO 

670 


IF (NX(l) .GT. NX<2) ) C0RE3= 0*0 


SLO 

680 


IF (NX(6) .LT. NX(5) ) C0RE2= 0.0 


SLO 

690 


IF (CORE3 .EO, 0,0) GO TO 102 


SLO 

700 

C 



SLO 

710 

C 

tertiary CORE REGION 


SLO 

720 

C 



SLO 

730 


NPRR=NX(U*l 


SLO 

740 


DO 70 I=MPRR*N 


SLO 

750 


IF (DIP (I) -GT. T0L2) go to 72 


SLO 

760 

70 

CONTINUE 


SLO 

770 

72 

CONTINUE 


SLO 

780 


NX(2»- I-l 


SLO 

790 


NPRRsNX (2) ♦! 


SLO 

800 


DO 75 I-NPRR*N 


SLO 

810 


IF (OIP(I) .LT, TOL2) GO TO 77 


SLO 

820 

75 

CONTINUE 


SLO 

830 

77 

CONTINUE 


SLO 

840 


NX<3)s I-l 


SLO 

850 

102 

CONTINUE 


SLO 

860 


IF (CORE2 .EO. 0,0) GO To 101 


SLO 

870 


93 



c 

C SECONDARY CORE REGION 

C 

MPRR-NX (6) 

r>0 80 M=?tNPRR 

1= NX(6) ♦ (?-M) 

IF (DlPin .GT. TOL3) GO TO 82 
80 CONTINUE 

82 CONTINUE 

MX(5)= I 
NPRR=NX <5) 
no 85 M=?,NPRR 
1= NX (5) ♦ (2-M) 

IF (DlPtl) .LT. T0L3) GO TO 88 
85 CONTINUE 

88 CONTINUE 

MX(4)= I 
IGl CONTINUE 
C 

C CHECK IF PRIMARY CORE EXISTS 

C 

IF (CORE! .EO. 0,0) GO TO 110 
IF (COPES ,NE. 0.0) GO TO 120 
no 130 I=2*N 

IF (DIP(I) ,LT. TOL2) GO TO 132 
130 CONTINUE 

132 CONTINUE 

NX(3)= I-l 
120 CONTINUE 

IF (CORE2 ,NE. 0.0) 60 To 140 

00 135 Ms2»n 
1= N ♦ (?-M) 

IF (DIP(I) .LT, TOL3) GO TO 137 
135 CONTINUE 

137 CONTINUE 

MX(4>= I 
140 CONTINUE 

IF (NX(3) ,GE. NX(4) ) CORE1=0.0 

110 CONTINUE 

C 

C CHECK VALUES OF * COREl* CORE2» CORE3 

C NhAX is MAXIMUM U{ ) VALUE 

C 

IF (COREl ,NE. 0.0 ) GO To 670 
C 

C COREl = 0.0 

C 

NX (3) =NMAX 
NX (4) =NMAX 

670 IF (C0RE3 .NE. 0.0> GO TO 680 
C 

C CORE3 = 0.0 

C 

NX (1 ) =NX (3) 

NX (2)=NX(3) 

680 IF (CORE? .NE. 0.0) GO TO 90 

C 

c Core? = o.o 

c 

NX (5)=NX(4) 

NX (6)=NX (4) 

90 CONTINUE 

NX1=NX ( 1 > 


5L0 880 
5L0 890 
SLO 900 
$LO 910 
JLO 920 
$L0 930 
fLO 940 
$LO. 950 
$LO 960 
SLO 970 
fLO 980 
$LO 990 
fLOlOOO 
fLOlOlO 
$L01020 
SL01030 
SL01040 
SL01050 
SL01060 
SLO1070 
SL01080 
5L01090 

sLonoo 

fLOl 1 10 

SL01120 

$L01130 

SLOH40 

fLOllSO 

SL01160 

fLOl 170 

fLOlieO 

SLOl 190 

SL01200 

fLO12l0 

SLO1220 

5L01230 

JLO1240 

5L01250 

$L01260 

$L01270 

SLO1280 

$L01290 

$L01300 

SL01310 

5LO1320 

SL01330 

ILO1340 

$L01350 

tL0l360 

tLOl370 

5L01380 

SL01390 

fL01400 

5L01410 

SL01420 

$LO1430 

fLO1440 

tL01450 

fL01460 

fL01470 

$LO1480 

SL01490 


94 



c 

c 

c 

c 

c 

c 

c 

c 

c 


c 

c 


2 


50 


52 


C 

C 

C 


NX2=NX (2) 

N*f3=NX (3) 

NX4=NX (4) 

NX5=NX (5) 

NX6=NX (6) 

C0REI=DN(NX4)-DN(MX3> 

CORE2rDN(NX6)*-DN'(NX5) 

COPE3=ON(MX2)-DN{MX1) 

0FL1=0N(NX1 ) 

DFL3«0N{NX5)-DN(NX4) 

DEL2=nN(NX3)-DN(Nx2> 

DEL4=DN<N)-0N(NX6> 

IF ( (COREl.EO.O.) . and. (CORE2.EO.O.) .AND. (CORE3.EO.O.) ) NCORE^ 0 

RP-TURN 

End 

subroutine cords (N.X*Y*THETA*U.RH0.PSI»DN»N0UCT,D0IST) 


subroutine cords 

CALCULATES 

DISTANCE BETWEEN STREAM LINES (ARRAY DN) 

DRAWS ARCS FROM TOP AND BOTTOM WALL TO MIDDLE OF DUCT 
compares Y COORD OF POTH ARCS AT MIDDLE POINT 

«««»*•«••«•»«»«»«•«*«»««»»«*»«««««»««•'«««* 

RFA'L M(190) 

DIMENSION X(190) *Y(190) ♦ PO ( 1 90 1 » H ( 190 I ♦ THETA (1 90 » 2) »R(190»2) * 

1 U(190*2) *TSn90,?) »RHO(190*2) 

DIMENSION OMG(190> *DN(190) 

DOUBLE OMG(190) *ON(190> 

DIMENSION PSI (IBO) 

DOUBLE PSI(190> 

ONd) * 0. 

DO 2 J=2*N 
JM=J-1 

7=2./(RH0( JMt2)*U(JM,2)*RHO(J»2)*U( J»2) ) 

ON(J)= ON(JM) ♦ Z * (PSI(J) - PSl(JM>) 

CONTINUE 

DO 50 Js=2»NDUCT 
JM=J-1 

X (J)=: X< JM)-(DN<J)-DN(JM) )*S IN (THETA (J»2) ) 

Y( J)= Y ( JM) ♦ (DN( J)-DN(JM) )*C0S(THETA( J,2) ) 

CONTINUE 

MPARsNDUCT^I 

KipAPR=N-l 

no 52 MM=NPAR,*NPARR 

jsN-MM*NDUCT 

JMssJ-il 

X (J>= X (JM)-(DN( J)-DN( JM) > «SIN( THETA (J»2> > 

V (J)= Y(JM)*(DN(J)-DN(JM) )*COS (THETA (J*2> ) 

CONTINUE 

YCOMP= Y(NDUCT*l> ♦ (DN(NDUCT) - DN(NDUCT*1>) * COS ( THET A ( J » 2 ) ) 
DDIST= (Y(NDUCT) - YCOMP) / DN(N) 
return 
End 

subroutine remove (N»NN*NGR lD*NX.NMAX»OM6t SI *S2*S3f S4»S5»DS*VIS* 

I ON *X*Y.M*P0*H, THETA* R,PS»U*TS*RHO) 

SUBROUTINE REMOVEpTS 


5L01500 

$L01510 

SL01520 

5L01530 

$L01540 

SL01550 

SL01560 

SL0I570 

SL01580 

SLO1590 

SL01600 

SL01610 

JL01620 

5LO1630 

SL01640 

SL01650 

SCO 

0 

SCO 

10 

SCO 

20 

SCO 

30 

SCO 

40 

SCO 

50 

SCO 

60 

SCO 

70 

SCO 

80 

SCO 

90 

SCO 

100 

SCO 

no 

SCO 

120 

SCO 

130 

SCO 

140 

SCO 

150 

SCO 

160 

SCO 

170 

SCO 

180 

SCO 

190 

SCO 

200 

SCO 

210 

SCO 

220 

SCO 

230 

SCO 

240 

SCO 

250 

SCO 

260 

SCO 

270 

SCO 

280 

SCO 

290 

SCO 

300 

SCO 

310 

SCO 

320 

SCO 

330 

SCO 

340 

SCO 

350 

SCO 

360 

SCO 

370 

SCO 

380 

SCO 

390 

SRE 

0 

SRE 

10 

SRE 

20 

SRE 

30 

SRE 

40 


95 


O U o (J.U O ID 


c 

c 

c 

c 

C 

C 

C 

C 

C 

C 

C 

C 

C 

C 

C 


c 


300 


1 
50 





SRE 

50 

REMOVES GPIO POINTS AS PROGRAM ITERATES 



5RE 

60 

(TO REDUCE COMPUTING TIME) 



5RE 

70 

REQUIRED parameters: 



?RE 

80 

N* NGRID* ARRAY U(90) 



$RE 

90 

(all other ARRAYS) 



SRE 

100 

method: 



SRE 

110 

FOR EVERY TENTH NSTEP 



$RE 

120 

SCANS ARRAY U 



IRE 

130 

IF (U(I-l) - U(I) / U(I) ) LESS THAN TOLERANCE* 


IRE 

140 

THEN GRID POINT IS REMOVED 



SRE 

ISO 

CORE REGIONS ARE NOT SCANNED 



SRE 

160 

arrays copied ape (0MG*S1 *S2*S3*S4*S5*DS*VIS*DN*X,Y*M,P0*H 

9 

SRE 

170 

THETA,R*PS*U»TS*RH0) 



SRE 

160 




SRE 

190 




SRE 

200 

REAL M(190) 



SRE 

210 

Dimension si (19o) *s?(19o) »s3U90) *54(190) *s5(190)*e<190) » 

El (190) * 

SRE 

220 

1 OS(190) «VIS(190) ,X(190) * Y ( 190 ) *P0 ( 190 ) *H ( 190 ) * 



SRE 

230 

2 THETA (190 *2) *R(190*2) *U ( 190 *2) * TS ( 190*2) *RHO(190*2) * 

SRE 

240 

3 NXSAVE(6) *NX(6) 



SRE 

250 

DIMENSION OM6(190) *ON( 190), PS (190*2) 



SRE 

260 

DOUBLE OHG(190) *ON(190) *PS(190*2) 



SRE 

270 

NGRIDs NGRI0*1 



SRE 

280 

IF (NGRID .LT* 10) RETURN 



SRE 

290 

NGPIC= 0 



SRE 

300 

DO 300 1=1*6 



SRE 

310 

NXSAVE(I>= NX (I) 



SRE 

320 

CONT INUE 



SRE 

330 

Js 1 



SRE 

340 

Do 24 I=2«N 



SRE 

350 




SRE 

360 

DO NOT REMOVE GRID POINTS FROM CORE REGIONS 



SRE 

370 




SRE 

380 

IF MI.GT. (NX5AV£(1)-3) ) .AND. ( I .LE. (NXSAVE (2) *3) ) ) GO 

TO 

51 

SRE 

390 

IF (d.GT. (NXSAVE(3)-3) ) .ANO. ( I .LE. (NXSAVE (4) +3) ) ) 60 

TO 

51 

SRE 

400 

IF (d.GT. (NXSAVE(5)-3> ) .AND. ( I ,LE. (NXSAVE (6) *3) ) > GO 

TO 

51 

SRE 

410 

IF (I .EQ, N) 60 TO 51 



SRE 

420 

URs ABS((U(I*2) • U(J*2>) / U(I*2)) 



SRE 

430 

ONN= ABS( (ON(I) - ON(J) ) / 0N(N) ) 



SRE 

440 

IF ( (UR. .LT. 0.02) .AND. (DNN .LT. 0.015)) GO TO 5*0 



SRE 

450 




SRE 

460 

GO TO STATEMENT SO WILL REMOVE GRID POINT 



SRE 

470 




SRE 

480 

J= J ♦ 1 



SRE 

490 

OMG(J) = OMG(I) 



SRE 

500 

DS(J) “ DSd) 



SRE 

510 

VTS(J) ViSd) 



SRE 

520 

On(J) = ONd) 



SRE 

530 

X(J) = X(I) 



SRE 

540 

Y(J) = Yd) 



SRE 

550 

M(J) = Md> 



SRE 

560 

P0(J) = POd) 



SRE 

570 

H(J) = Hd) 



SRE 

580 

ThETA(J.2) = THETA(I*2) 



SRE 

590 

R(J*2) = Rd*2) 



SRE 

600 

PS(J.?> = PSd,2) 



SRE 

610 

U(J*2) = Ud*2) 



SRE 

620 

TS(J*2) = TSd.P) 



SRE 

630 

RhO(J* 2) = RHOd*?) 



SRE 

640 

IF (NXd) .EO. I) NX(1> s J 



SRE 

650 

IF (NX(2) .EO. I) NX(2) = J 



SRE 

660 

It, <NX(3) .EQ. I) NX(3) = J 



SRE 

670 
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fiO TO 320 

*AH 

450 

350 

CONTINUE 

SAR 

460 


HH/k-XK 

SAR 

470 


nRCA=ORC 

SAP 

480 


IF(FW*DRC .GT.0.0) DRC8* DRCR*0.5 

SAR 

490 

320 

CONTINUE 

SAR 

500 


FW=* DRC 

SAR 

510 

C 


SAR 

520 

C 

CHECK FOR TERMINATION 

SAR 

530 

C 


SAR 

540 


IF(ABSIXXB-XXA) .LF. l.E-6) GO TO 310 

SAR 

550 


1F(ARS(0RC) ,LE* l.E-9 ) GO TO 310 

SAR 

560 

300 

CONTINUE 

SAR 

570 


WRITE<6*900) 

SAR 

580 

900 

FORMAT </7Xt3SHN0 CONVERGENCE IN FINDING XC AND YC/I 

SAR 

590 

310 

CONTINUE 

SAR 

600 


RETURN 

SAR 

610 


End 

SAR 

620 


PEAL FUNCTION CENTRE ( XNOZ *YN0Z*MN02 *BNOZ*XX »RC*XW*YW*NPAIR*J) 

SCE 

0 

C 


5CE 

10 

C 


SCE 

20 

C 

FUNCTION CENTRE 

SCE 

30 

C 

AIDS IN FINDING A SUITABLE TANGENT TO UPPER WALL 

SCE 

40 

C 

GIVEN XNOZ* XX - PROJECTS TWO TANGENTS 

SCE 

50 

C 

FINDS INTERSECTION XC, YC 

SCE 

60 

C 

returns nlFFFRENCE BETWEEN LENGTH OF RADIUS 

SCE 

70 

C 


SCE 

80 

c 


SCE 

90 


PEAL MN07* MM 

SCE 

100 


dimension XW(99*2) tYW(99*2) 

SCE 

110 


DIST (XX*YY*XXC»YYC)=SORT ( ( ABS < XX-XXC) **2> ♦ (ABS(YY-YYC) **2> ) 

SCE 

120 


CALL WALLS (XW*YW*XX*VY*TT,RR»NPAIR,J) 

SCE 

130 


TAN(TT) 

SCE 

140 


RB= yY - (MM* XX) 

SCE 

150 


IF (ABS (MNOZ-MM) .GT, O-OOOOl) XC= (80-BNOZ) / (MNOZ-MM) 

SCE 

160 


IF (ABS(MNOZ-MM) .LE. O.OOOOl) XC= (XNOZ^XX) * 0,5 

SCE 

170 


YC= mnoz*xc ♦ bnoz 

SCE 

180 


RC= DIST (XNOZ*YNOZfXC*YC) 

SCE 

190 


centre* OIST(XNOZ,YNOZ*XC*YC) - DiST (XX» YY»XC, YC) 

SCE 

200 


RETURN 

SCE 

210 


Evf) 

SCE 

220 


SUBROUTINE OMGSET (OMG, JSTaRT , JF IN, jSTSER f JFNSER ,DMST tDMFNf OMM lOt 

SOM 

O' 


1 RSTjRFIn) 

SOM 

ICi 

c 


SOM 

20 

c 


SOM 

30 

c 

SUBROUTINE OMGSET 

SOM 

40 

c 

SETS OMG DISTRIBUTION ARRAY OMG 

SOM 

50 

c 

FOR EACH CORE REGION, OMGS ARE CONSTANT IN MIDDLE 

SOM 

60 

c 

CLOSE together AT EDGES 

SOM 

70 

c 

OMG(JSTART) and 6MG(JFIN> must already be SPECIFIED 

SOM 

80 

c 


SOM 

90 

c 


SOM 

100 


DIMENSION OMG(190) 

SOM 

no 

c 

Double omgu9o> 

SOM 

120 


JC0NST= JSTaRT ♦ JSTSER ♦ 1 

SOM 

130 


riMG.( JCONST)* OKG(JSTART) ♦ DMMID 

SOM 

140 


NPRR=JC0NST*1 

SOM 

150 


NPRRR* JF IN- JFNSFR- 1 

SOM 

160 


Do 500 J*NPRR,NPRPR 

SOM 

170 

500 

OMG( J)sOMG< J-1) ♦DMMID 

SOM 

180 

c 

SOL%E FOR RST IN SERIES A+AR+ AR**^* , . . 

SOM 

190 


XA* 1,5 

SOM 

200 


FXA* SERIES (DMMlDtDMST* JSTSER, XA) 

SOM 

210 


97 




IF (NX(4) ,eO, I) NX(4) = J 

SRE 

680 


IF (NX(5) *EO. I) NX(5) s J 

SRE 

690 


IF (MX(6) .EO» I) NX(6) 3 J 

SRE 

700 


IF (NMAX ,EO. I) NMAX* J 

SRE 

710 

24 

Continue 

SRE 

720 


N« J 

SRE 

730 


NN s N - 1 

SRE 

740 


00 20 Js2fNN 

SRE 

750 


JP= J+1 

SRE 

760 


JM= J-l 

SRE 

770 


S1(J> s OMGIJP) - OMG(JM) 

SRE 

780 


S2(J»= OMG(JP) - OMGU) 

SRE 

790 


S3(J)= OMGU) - 0H6UH) 

SRE 

800 


S4(J)s S2(J) / S3(J) / SKJ) 

SRE 

810 


S5(J>» S3(J) / S2U) / Sl(J> 

SRE 

820 

20 

CONTINUE 

SRE 

830 


rfturn 

SRE 

840 


E^'D 

SRE 

850 


SUBROUTINE ARCOIS (XNOZ* VNOZ*TNOZ*XX*RC*XW*YW*NPaIR* J) 

SAR 

0 

C 


SAR 

10 

C 


SAR 

20 

c 

SUBROUTINE ARCOIS 

SAR 

30 

c 

FINDS XX <ON wall SPECIFIED) SUCH THAT 

SAR 

40 

c 

ARC RADIUS RC IS TANGfNT AT XX 

SAR 

50 

c 

AND PASSES THROUGH N07ZLE AT ANGLE TNOZ 

SAR 

60 

c 


SAR 

70 

c 


SAR 

80 

c 

FIND EQUATION OF LINE TANGENT TO XI 

SAR 

90 

c 


SAR 

100 


real mnoz 

SAR 

110 


DIMENSION XW(99,2) *YW(09*2) 

SAR 

120 


DIST(XX*YY*XXC*YYC)sSQRT( (ABS(XX-XXC>**2>+(ABS(YY-YYC)**2) ) 

SAR 

130 


MNOZ= tan (TNOZ) 

SAR 

140 


PNOZ= YNOZ - MNOZ * XNOZ 

SAR 

150 

c 


SAR 

160 

c 

CHOOPPING technique TO FIND CENTRE (X2A) .LT. 0.0 

SAR 

170 

c 

And centre(X2b) .gt. o.o 

SAR 

180 

c 


SAR 

190 


DO 200 K sl*NPAlR 

SAR 

200 


XXB= XW(KtJ) 

SAR 

210 


nRCB= CENTRE<XN02>YN0 Z*Mn 0Z*8N0Z, XXB»RCfXW*YW*NPAIR,J) 

SAR 

220 


IF (DRCB .GT. O.Cf) 60 TO 210 

SAR 

230 

900 

CONTINUE 

SAR 

240 

210 

CONTINUE 

SAR 

250 


XXA= XW(K-2fJ) 

SAR 

260 


DRCAs CENTRE (XNOZ* YNOZfMNOZfBNOZvXXAtRC^XWvYWfNPAIRfJ) 

SAR 

270 


FW= ORCA 

SAR 

280 

C 


SAR 

290 

C 

MODIFIED REGULA FALSI ALGORITHM 

SAR 

300 

c 


SAR 

310 


00 300 K=l*50 

SAR 

320 

c 


SAR 

330 

c 

GUESS NEW X2 

SAR 

340 

c 


SAR 

350 


XXs (ORCB*XXA - DRCA*XXR) / (ORCB-DRCA) 

SAR 

360 


DRC= CENTRE (XNOZ* YNOZ* mnoz *BNOZ* XX* RC*XW*YW*NPAIR*J) 

SAR 

370 

r» 

o 


SAR 

380 

r* 

u 

CHANGE APPROPRIATE ENDPOINT TO MIDPOINT 

SAR 

390 

k# 


SAR 

400 


IF (ORCA*DRC .GT.0.0) GO TO 350 

SAR 

410 


XXB= XX 

SAR 

420 


DRCB»s DRC 

SAR 

430 


IF (FW*DPC.GT. 0.0) 0RCA»0RCA*.5 

SAR 

440 
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XP= 1*3 

50M 

220 


FXB= SEPIES(0MMID*0MST* JSTSER*XB) 

SOM 

230 


no 505 Jsl,50 

SOM 

2A0 


XR=s (FXB*XA - FXA*XB)/(FX8-FXA) 

SOM 

250 


IF (XB •RT, 2.0) XR= 2.0 

SOM 

260 


FXB= SERIES (OMMlO*DMSTtJSTSER»XB) 

SOM 

270 

c 


SOM 

280 

c 

CHECK FOR termination 

SOM 

290 

c 


SOM 

300 


IF (ARS(FXB) .LE. 1*E-04) GO TO 510 

SOM 

310 

505 

CONTINUE 

SOM 

320 


WRITE(6»900) JSTART 

SOM 

330 

900 

FORMAT (/7X*22HNO SOLUTION JSTART *fI5/) 

SOM 

340 

510 

CONTINUE 

SOM 

350 


RST» XB 

SOM 

360 

C 


SOM 

370 

C 

SOLVE FOR RFIN IN SERIES 

SOM 

380 

C 


SOM 

390 


XA= 1.5 

SOM 

400 


FXAs SERIES (DMM ID tOMFN*JFNSER*XA> 

SOM 

410 


1.3 

SOM 

420 


FXB= SERIES ( DMM TO, DMFN.JFNSER.XB) 

SOM 

430 


no 525 J=1.50 

SOM 

440 


XR= (FXB*XA - FXA*XB) / (FXB-FXA) 

SOM 

450 


IF <XB .GT. 2.0) X0= 2.0 

SOM 

460 


FXB= SERIES (DMM TO iDHFN.JFNSERtXB) 

SOM 

470 

c 


SOM 

480 

c 

CHECK FOP termination 

SOM 

490 

c 


SOM 

500 


IF (ABS(FXB) .LE. l.E-OA ) GO TO 520 

SOM 

510 

525 

CONTINUE 

SOM 

520 


WRlTF(6t901) JFIN 

SOM 

530 

901 

FORMAT (/7X.20HNO SOLUTION JFIN =.I5/) 

SOM 

540 

520 

CONTINUE 

SOM 

550 


RFIN= XR 

SOM 

560 


no 550 J=1.JSTSER 

SOM 

570 

550 

DMG( JSTART*J)= OMG ( JSTART ♦J-1 > ♦ DMST * (RST**(J-1)> 

SOM 

580 


no 560 J=1»JFNSEP 

SOM 

590 

560 

0MG(JFIN-J)= OMG(JFIn-J^I) - OMFN ♦ (RF IN** ( J-1 ) ) 

SOM 

600 


RETURN 

SOM 

610 


EN!D 

SOM 

620 


PEAL FUNCTION SEP lES ( SUM » A *N » R )- 

SSE 

0 


IF (R .NE, 1.0) SERIES* <R** (N* 1 ) -1 . 0 > / (R-1.0) - (SUM/A) 

SSE 

10 


IF (R .EQ. 1.0) series* N*1 - (SUM/A) 

SSE 

20 


RETURN 

SSE 

30 


END 

SSE 

40 
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Table 1 


Mixing Section Dimensions 
(Inches) 



LOWER 

WALL 


UPPLR 

WALL 


J 

X 

Y 

CURV 

X 

Y 

CUKV 

1 

-6.0000 

-4.1028 

.11979 

-3.3270 

6* 02J 6 

-.03137 

2 

-5.9500 

-4. 0266 

.12608 

-3.3000 

5.7680 

-.04041 

3 

-5. 9000 

-3.9523 

.13277 

-3.2500 

5.3572 

-.06887 

4 

-'5.5000 

-3.4J24 

,17553 

, -3.1250 

4.7764 

-.11339 

5 

-5.00C0 

-2.9342 

,21607 

-3,OOGO 

4,3953 

-.13194 

6 

-4.5000 

-2.5697 

.22500 

-2. 7500 

3. 8304 

-.17358 

7 

=T» . o 7iro 

“=?.2 98S 

T2 24i^ 

-2.5000 ~ 

3,4iJll 

-.20405 

d 

-3.5300 

-2. 1047 

•18402 

-2.2500 

3 ■ 0 j 3 

-.17665 

9 

-3.C00C 

-1.9628 

.15850 

-2. COCO 

2.7705 

-.16650 

10 

-2.5Qff0 

-1. 8637 

.13904 

^175000 

2.2963 

-. 18696 

11 

-2.3500 

-1. 8403 

,14043 

-l.OOCO 

1.9315 

-.20381 

12 

-2.3125 

-1.8349 

.14075 

-.5000 

1.6544 

-.20305 


-2.3125 

-1.8390 

.OOOOQ 

. C 000 

1 . 4449 

-.21‘235“ 

14 

-2.0000 

-1.7783 

.00000 

, 5000 

1.2971 

-.21607 

15 

-l.OOOC 

-1. 5840 

.0 00 00 

. 6250 

1.2689 

-.22019 

16 

.OOOO 

-1. 3897 

.0 oOdo 

. 7500 

1.2444 

-.22389 

17 

.7503 

-1.2440 

.00000 

. 7500 

1.2443 

-.08962 

18 

- - 

.7500 

-1.2443 

.08962 

. 80L0 

1.2350 

-.08984 


. S Q 0 0 

-1. 23515 

.15 8984 

.8500 

1.2260 

-.6915 tr7 

20 

.8500 

-1. 2260 

.09007 

1.0000 

1.2009 

-.07313 

21 

1.0000 

-1.2009 

,07313 

1.5000 

1.1298 

-.07108 

22 

1.5000 

-1.1298 

.D'71'08 

2.0000 

1 . 077'3 

-. 06664 

23 

2. GOOD 

-1. 0773 

.06664 

2. 5000 

1.0423 

-.05432 

24 

2.5000 

-1.0423 

.05432 

3.0000 

1.0205 

-.04063 

25 

3 • Q Q 0 G 

-1.13^05 

TO 4 a 63 

j;4OC50 

1 . 0 D^'5 

-.01971“ 

26 

3 .4000 

-1. 0095 

.0 1971 

3.4500 

1.0032 

-.01971 

27 

3.4500 

-1.0082 

.0 1971 

3.5000 

1.0 369 

-.01971 

~2B' 

3 . 5 0 0 0 

-1 . D1TB9 

• 0 19 71 

3 . 5 0 C 0 

r.Gu70 

. GOCOO'~ 

29 

3 .5000 

-1. 0070 

.00000 

4.0000 

. 9984 

.00000 

30 

4.0000 

-.9984 

.00000 

5. 0000 

. 9811 

.00000 

“31 

5^Tia(itr 

-."9<rrr ~ 

TO 0 0 013 

6.0000 

. 96T9 

706'0'09 

32 

6.0000 

-.9639 

.00000 

7.0000 

.9466 

. 00000 

33 

7, coco 

-.9466 

.00000 

7 . 5 C 0 0 

.9380 

. GQOOO 

34 

7.50'33 

-.9380 

TO joins ■ 

7750015 ■ 

.“938^0 

. C 0 0 0 0 

35 

7.50G0 

-.9380 

.OOODO 

6. 0000 

,9380 

.00000 

36 

8 .0 GOO 

-.9380 

.03000 

9.0000 

. 9380 

.00000 

3 f 

9.13 OTTO 

=7^310 

721300D 

rf.-coco 

,'93315 

7ca'oco“ 

38 

1C .0000 

-.9380 

.03000 

10.5000 

. 9380 

.00000 

39 

10 .5000 

-.9380 

.08000 

10.5000 

. 9380 

.00000 

40 

icr;5T(n; 

-.938H 

TCTO 0*0 

I'i.OOOO 

.■9642 “ 

.OCUQ'O' 

41 

11 .oooc 

-.9642 

.QGQOO 

12.0000 

1.0166 

.00000 

42 

12.0000 

-1. 0166 

.00000 

14.0000 

1.1214 

.00000 

43 

14. 0300 

“r.T?I4 

rOGODO 

16.0000 

~T72Z^2 

. 6 OlilS 0 

44 

16.0000 

-1.2262 

.00030 

16.0000 

1.3310 

.00000 

45 

18.0000 

-1.3310 

.00000 

20.0000 

1.4358 

.00000 

“" 4‘6 ' 

20.01T3 0 

-1. 4 3 58 

. 00000 ' 

2 2 • 0 C 0 0 

1. 5406 

• DTCQU 

47 

22.0030 

-1. 5406 

.00003 

23. COCO 

1. 5930 

*00000 
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Table 2 

Variation of Individual Integrated Traverse Mass 
Flows for Each Test Run 
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Table 3 


Traverse Locations and Data Summary 


Maximum Flow - Atmospheric Discharge 


Reduced Flow - Back Pressure at Discharge 


Test Number 

9 1 Q 

! 2 3 

6 

z 

1 

4 

Ut 

00 

Nozzle Position (See Fig. 3) 

r 1 

2 

3 

4 ! 1 

2 

3 ^ 4 

Nozzle Angle 

22.5° 

45° 

45° 

1 

67.5° i 22.5° 

45° 

45° 

67.5° 

Nozzle meter 

Spacing inch 


0.020 

0.80 

0.034 

1.32 

0.018 

0.70 

0.020 

0.80 

o 

<N O 
O 00 

o o 

0.034 

1.32 

0.018 

0.70 

Traverse Location + 









s = -0.070 meters, - 2.75 inches 




P 

max 




P 

max 

s = -0.024 meters, - 0.95 inches 


P 

max 

P 

max 

P, P 

max 


p 

max 

P 

max 

P, P 

max 

X = +0.013 meters, + 0.50 inches 

P 1 P 

max ^ max 

P 

max 

P 

max 

P 

max 

P 

max 

P 

max 

P 

max 

X = +0.038 meters, +1.50 inches 





P 

max 

P 

max 

P 

max 


X = +0.064 meters, + 2.50 inches 









X = +0.114 meters, +4.50 inches 

U, P 

max 

U, P 

max 

U, P 

max 

U, P 

max 

U, P 

max 

U, P 

max 

U, P 

max 

U, P 

max 

X = +0.165 meters, +6.50 inches 




■BiiH 

■Bini 

U, P 

max 

U, P 

max 

U, P 

max 

X = +0.254 meters, +10.00 inches 

U, P 

max 








X = +0.318 meters, +12.50 inches 


U, P 

max 

U, P 

max 

U, P 

max 

U, P 

max 

U, P 

max 

U, P 

max 

U, P 

max 

X = +0.419 meters, +16.50 inches 




U, P 

max 



l_ 


X = +0.521 meters, +20.50 inches 

U, P 

max 

U, P 

max 

U, P 

max 

U, P 

max 

U, P 

max 

U, P 

max 

U, P 

max 

U, P 

max 


Axial Wall Static Pressure Profiles (Figures 16 to 19) 
P = Vertical Total Pressure Profile (Figures 21 and 22) 
^max ~ Maximum Pressure Location (Figures 23 to 30) 

U = Vertical Velocity Profile (Figures 31 to 38) 
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Table 4 


Summary of Experimental Test Conditions and Flow Rates 

2 

Nozzle Throat Area = .9688 in 

Mixing Section Throat Size = 1.875 inch 

Nozzle Pressure = 36.70 psia (constant) 














































Table 5 


Comparison of Experimental and Analytical Mass Flow Rates 


Run 

No. 

Mixing 
Section 
Mass Flow 
Rate From 
Traverse Data 
lb/ (s^-m) 

HiiliSilil 



Comparison 

of 

Traverse to 
Analytical 
Mass Flow 

^ir 

Comparison 

of 

Orifice to 
Analytical 
Mass Flow 

@“(I) 

(2) 

1 

.4995 

_ 

— 

.4905 

+ 1,8% 


2 

.4332 

.3776 

+ 12.8% 

.4086 

+ 5.7% 

- 8.2% 

3 

.4627 

.3954 

+ 14.5% 

.4175 

+ 9.8% 

- 5.6% 

4 

.5238 

- 

- 

.4930 

+ 5.9% 

- 

5 

.5208 

- 

- 

.5100 

+ 2.1% 

- 

6 

.4220 

.3886 

+ 7,9% 

.4117 

+ 2.4% 

- 5.9% 

.■ -7 

.4390 

.3846 

+ 12.4% 

.3960 

+9.8% ■ 

- 3.0% 

8 

.5242 

- 

- 

.4840 

+ 7.7% 

- 


O 

Ln 









































Primary Nozzle 
Side Plates 

Nozzle Positioning Plates 
Suction Plenum 
Upper Contoured Plate 
Lower Contoured Plate 
Plexiglass Windows 
Inlet Bellmouth 
Coanda Surface 


Figure 2 

Assembly Sketch of Two Dimensional Ejector Test Rig 
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Figure 4 
Primary Nozzle 
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1. Nozzle 

2. Mixing Section Side Plate 

3. Top Contoured Plate 

4. Bottom Contoured Plate 

5. Table Top 

6. Screen 




Figure 6 

Side View of Test Section Showing 
Nozzle Positioned at 

22.5® and 0.020 meters (0.80 inches) spacing 
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Figure 7 
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Side View of Test Section Showing 
Nozzle Positioned at 

45® and 0.034 meters (1.32 inches) spacing 



Figure 8 

Side View of Test Section Showing 
Nozzle Positioned at 

67.5® and 0.018 meters (0.7 inches) spacing 
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Boundary Orifice for 


Layer Suction suction flow 
System 


Figure 9 


Schematic of Experimental Layout 
for Reduced Flow Conditions 









Figure 10 

View of Test Facility Showing Plenum, 
Mixing Section, Manometer Board, and 
Pressure Sampling Valves 
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Figure 11 

Top View of Mixing Section Showing 
Static and Traversing Tap Locations 
Traversing Probe is Positioned at 
"A" Location on Coanda Surface 
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Figure 14 

Side Viex-j of Test Section Shovjing 
Extended Inlet, lytLxing Section and Suction 
System Piping and Manifold 
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Wall Static Pressure - inches of water gage (x = 10.0 inch) 


Mixing Section Total Mass Flow = kg/sec. -m. 
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Figure 15 

Comparison of Experimental and Analytical Flow Rates 


Wall Static Pressure - kPa gage 






Wall Static Pressure - Inches of Water Gage 
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Solid Lines are Analytical Results 




Figure 16 


Comparison of wall static pressure distributions 
in a s}iTnmetrical mixing section with a 1.875” throat 
computed with the present streamline coordinate pro- 
gram and with the program from CR-2251 
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Wall Static Pressure - Inches of water gage 


Run 

Data Analytical 
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A Top 


# Bottom 
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A Top 


O Bottom 
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Figure 17 Wall static pressure distributions for Tests 1 and 2. 

Nozzle positioned at 22.5® and 0 , 020 meters (0.80 inches) 
spacing 


Wall Static Pressure - kPa gage 




Wall Static Pressure - inches of water gage 





Nozzle positioned at 45° and 0.020 meters (0.80 inches) 
spacing 
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Wall Static Pressure - kPa gage 


Wall Static Pressure - Inches of water gage 



Distance Along Distance From Leading Edge of Contoured 

Curved Surface Plates 



Nozzle positioned at 45® and 0.034 meters (1.32 inches) 
spacing 
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Wall Static Pressure - kPa gage 


Wall Static Pressure - Inches of water gage 
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a-^ — f— ►x meter 

-0.10 0 0.10 0.20 0.30 O.AO 0.50 0.60 



Curved Surface Plates 




Nozzle positioned at 67.5® and 0.018 meters (0.70 inches) 
spacing 
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Wall Static Pressure - kPa gage 


Wall Static Pressure - inches of water gage 



VJall Static Pressure - kPa gage 









Figure 23 

Locii of Maximum Stagnation Pressures (P ) in the Mixing Section 

° o max 

Nozzle Positioned at 22.5® and 0.020 meters (0.80 inches) spacing 
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Figure 24 

Lodi of Maximum Stagnation Pressure (P ) in the Mixing Section 

o max 

Nozzle Positioned at 22.5“ and 0.020 meters (0.80 Inches) spacing 
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Figure 25 


Locii of Maximum Stagnation Pressure (P ) in the Mixing Section 

® o max ® 

Nozzle Positioned at 45® and 0.020 meters (0.80 inches) spacing 
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Figure 27 


Locii of Maximum Stagnation Pressure (P ) in the Mixing Section 

o max 

Nozzle Positioned at 45° and 0,034 meters (1.32 inches) spacing 
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Test 6 

O Experimental 

Analytical (Starts at the nozzle) 
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Figure 28 


Locii of Maximum Stagnation Pressure (P ) in Mixing Section 

o max 

Nozzle positioned at 45® and 0.034 meters (1.32 inches) spacing 
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Figure 29 

Locii of Maximum Stagnation Pressure (P^ Mixing Section 

Nozzle positioned at 67.5° and 0,018 meters (0.70 inches) spacing 
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Figure 31 Total Pressure Profiles for Run 7 at x = +0.013 

meters (+0.50 inches). Nozzle Positioned at 67.5®, 
0.018 meters (0.70 inches) spacing 







o max X = constant 
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Figure 32 Total Pressure Profile for Run 8 at x = +0,013 

meters (+0.50 inches). Nozzle Positioned at 67.5®, 
0.018- meters (0.70 inches) spacing 
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c = + 0.520 meters (+ 20.50 inches) 


X = + 0.420 meters (+ 16.50 inches) 


X = + 0.318 meters (+ 12.50 inches) 
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Figure 33 

Velocity Profiles for Run 1 , 

Nozzle 

Positioned at 


22.5® and 0.020 meters (0.80 

inches) 

spacing 







wall 


wall 


Figure 34 Velocity Profiles for Run 2 ^, Nozzle Positioned at 
22.5° and 0,020 meters (0.80 Inches) spacing 
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max X = constant 
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constant 



Figure 36 Velocity Profiles for Run Nozzle Positioned at 
45° and 0.020 meters (0.80 inches) spacing 
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IVF I 





max X = constant 



Figure 37 Velocity Profiles for Run _5, Nozzle Positioned at 
45° and 0.034 meters (1.32 inches) spacing 
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maxjx = constant 
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Figure 38 Velocity Profiles for Run Nozzle Positioned at 
45° and 0.034 meters (1,32 inches) spacing 
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ma^yx = constant 
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Figure 39 Velocity Profiles for Run Nozzle Positioned at 
67.5® and 0.018 meters (0.70 inches) spacing 




max X = constant 
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Figure 40 Velocity Profiles for Run Nozzle Positioned at 
67.5° and 0.018 meters (0.70 inches) spacing 
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